{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "97Ad8yROpHV8"
      },
      "source": [
        "# Filtering Signals\n",
        "In this tutorial we will practice filtering MIMIC waveform signals.\n",
        "\n",
        "Our **objectives** are to:\n",
        "- Filter signals using the [SciPy signal processing package](https://docs.scipy.org/doc/scipy/reference/signal.html).\n",
        "- Understand how to interpret the amplitude-response of a filter.\n",
        "- Gain experience in filtering PPG signals."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "efztffyOpHV-"
      },
      "source": [
        "<div class=\"alert alert-block alert-warning\"><p><b>Context:</b> Filtering is used to eliminate noise from physiological signals. For instance, ECG signals can contain mains frequency noise due to electrical interference. Ideally, a filter would attenuate unwanted frequency content in a signal whilst retaining the physiological frequency content.</p></div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1y92CrBlpHV_"
      },
      "source": [
        "<div class=\"alert alert-block alert-info\"><p><b>Extension:</b> If you've not seen it before, then have a look at the <a href= https://docs.scipy.org/doc/scipy/reference/signal.html>SciPy signal processing package</a>. How might it be helpful for processing PPG signals?</p></div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "McluxdGrpHV_"
      },
      "source": [
        "## Setup"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "742rlU3ApHV_"
      },
      "source": [
        "_The following steps have been covered in previous tutorials. We'll just re-use the previous code here._"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "byobDIEIpHV_"
      },
      "outputs": [],
      "source": [
        "# Packages\n",
        "import sys\n",
        "from pathlib import Path\n",
        "!pip install wfdb==4.0.0\n",
        "import wfdb\n",
        "\n",
        "# The name of the MIMIC IV Waveform Database on Physionet\n",
        "database_name = 'mimic4wdb/0.1.0'\n",
        "\n",
        "# Segment for analysis\n",
        "segment_names = ['83404654_0005',\n",
        "                 '82924339_0007']\n",
        "\n",
        "segment_dirs = ['mimic4wdb/0.1.0/waves/p100/p10020306/83404654',\n",
        "                'mimic4wdb/0.1.0/waves/p101/p10126957/82924339']\n",
        "\n",
        "rel_segment_n = 0\n",
        "rel_segment_name = segment_names[rel_segment_n]\n",
        "rel_segment_dir = segment_dirs[rel_segment_n]\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iUNUDMnopHWB"
      },
      "source": [
        "---\n",
        "## Extract one minute of PPG signal from this segment"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A0YsdFikpHWB"
      },
      "source": [
        "_These steps have been covered in previous tutorials, so we'll just re-use the code here._"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "2geoFoDBpHWC",
        "outputId": "29847fb7-2f05-4762-b311-e40ce5a6136c",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Metadata loaded from segment: 83404654_0005\n",
            "60 seconds of data extracted from: 83404654_0005\n",
            "Extracted the PPG signal from column 4 of the matrix of waveform data.\n"
          ]
        }
      ],
      "source": [
        "# time since the start of the segment at which to begin extracting data\n",
        "start_seconds = 20 \n",
        "n_seconds_to_load = 60\n",
        "\n",
        "segment_metadata = wfdb.rdheader(record_name=rel_segment_name, pn_dir=rel_segment_dir) \n",
        "print(f\"Metadata loaded from segment: {rel_segment_name}\")\n",
        "\n",
        "fs = round(segment_metadata.fs)\n",
        "sampfrom = fs*start_seconds\n",
        "sampto = fs*(start_seconds + n_seconds_to_load)\n",
        "\n",
        "segment_data = wfdb.rdrecord(record_name=rel_segment_name,\n",
        "                             sampfrom=sampfrom,\n",
        "                             sampto=sampto,\n",
        "                             pn_dir=rel_segment_dir)\n",
        "\n",
        "print(f\"{n_seconds_to_load} seconds of data extracted from: {rel_segment_name}\")\n",
        "\n",
        "for sig_no in range(0, len(segment_data.sig_name)):\n",
        "    if \"Pleth\" in segment_data.sig_name[sig_no]:\n",
        "        break\n",
        "\n",
        "ppg = segment_data.p_signal[:,sig_no]\n",
        "fs = segment_data.fs\n",
        "print(f\"Extracted the PPG signal from column {sig_no} of the matrix of waveform data.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "beaUjk5spHWC"
      },
      "source": [
        "---\n",
        "## Create a filter"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "scLn9bVKpHWD"
      },
      "source": [
        "- Import the [SciPy signal processing package](https://docs.scipy.org/doc/scipy/tutorial/signal.html), which contains functions for filtering and differentiating."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "1db-InU9pHWD"
      },
      "outputs": [],
      "source": [
        "import scipy.signal as sp"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MuzX0XoCpHWD"
      },
      "source": [
        "- Specify the high- and low-pass filter cut-offs"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "cGKF92QTpHWD"
      },
      "outputs": [],
      "source": [
        "# Specify cutoff in Hertz\n",
        "lpf_cutoff = 0.7 \n",
        "hpf_cutoff = 10"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RtLUEaAHpHWE"
      },
      "source": [
        "- Create a Butterworth filter using the [butter](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html) function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "HmiOvyHZpHWE"
      },
      "outputs": [],
      "source": [
        "sos_ppg = sp.butter(10,\n",
        "                    [lpf_cutoff, hpf_cutoff],\n",
        "                    btype = 'bp',\n",
        "                    analog = False,\n",
        "                    output = 'sos',\n",
        "                    fs = segment_data.fs)\n",
        "\n",
        "w, h = sp.sosfreqz(sos_ppg,\n",
        "                   2000,\n",
        "                   fs = fs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "k84fKyK6pHWE"
      },
      "source": [
        "- Plot filter characteristics"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "id": "fjd8i8-SpHWE",
        "outputId": "951c910b-e786-4598-c568-6064effb9d71",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        }
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3gc1dXA4d+RrGbJltw77ja44CbbQCDYdAjBhFADuILpEEJCCfkIJLSEAAk9Nu40O/QeqiEUdwzYGFwB495tWV063x9zZdZCWq2s3Z1d7XmfZx/tzszOnJ0d3bP3zp07oqoYY4wx4ZbkdwDGGGPqJ0swxhhjIsISjDHGmIiwBGOMMSYiLMEYY4yJCEswxhhjIsISjDlgIjJMRH6oxfKzReSiSMbkx7ZqIiIZIvKKiOwSkf+IyPki8lbAfBWRbtHafqS2Y0xllmDqQES+FZECEckTkR0i8pqIdKjF+/crWGpbYEdbpAvCeuxMoBXQTFXPUtUnVfWEqhYUkakicnsktx/mdRtTLUswdfdLVc0C2gCbgAf9CkREGsTTehNIR2C5qpZGekPVfFdBt59o32+ifV5fqao9DvABfAscF/D6FLx/5IrXs4GLAl6PBj5yzz8EFNgL5AGjgAKg3L3OA9ri/Qi4EVgFbANmAU3dOjq5dYwDvnfrnAZc5+a3c/OvcK+7AtuBJPf6YmClm/Yy0DYgVgWuAFYAa6qI9xxgGPADcB2wGdgAjAmyv2YDdwHzgN3ASxWfxc3/D7AR2OW21ztg3lTgYeA1YA8wF+gaMP944Gv33oeADyr2vdvvH7vpu9xyxwa8dwywzK13NXBJwLzmwKvATref/hew/24A1rn3fRO4zoD33wYUAyVuv40LPA4C9nU3YLxbrtgt+4qb3xZ4DtjivourA957K/As8ITbpxeFuP2PgfvxjqnbgTTgH3jH0SbgMSAjYD1/cN/vemBsRcw1Hefu9cHA227/fQOcXYvvtXfAezcBfwRaA/l4NbKK5Qa6/ZNSxXfwk30EZAOT3Gda5/ZBslu+G97xswvYCsys9F1d7Y6TrcA9AcdDEvAn4Du8/4fpQHal/9VRbh9vBW4OWO8QYIGLbxNwX8C8w4BP8I7Bz4Fhfpd9IZeRfgcQzw8CEgzQEK9wnx4wv6Z/vH3/pO71MOCHStu4BpgDtHeFwL+Bp928ioN2OpAJZLh//oqC6Td4iWmmez0WeMk9P8Yd5APdeh8EPqwU29tAU1xBU028pcBfgBS8BJsPNKlmf812/8x9XLzPAU8EzB8LNHLx/BNYHDBvKl5hOARoADwJPOPmNccrnM50cVzr4gpMMKVuegpectzFj4n6F3jJV4Cj3WcY6ObdhVfYprjHUW65nsBaXFJ230XXaj73rZU+Z7XHgfuctwfMSwIWArcAqUAXvMLtxIB1lwCnu2UzQtx+KXCV25cZeMnmZfd9NwJeAe5yy5+EV+hVfG9PEWKCccuvxUviDYABeMddrxC+10Z4CeA6IN29HurmvQ5cFrDN+4EHg+z//fYR8ALe/1Im0BLvR88lbvmngZvdsunAkZW+q/fdfjoIWM6Px9lYvB9sXYAs4HlgRqX/1Ylu+/2AIuAQN/9T4EL3PAs4zD1v5/bPKS6e493rFn6XfyGVkX4HEM8PvASTh/fLogTv113fgPnV/uO516EkmGXs/2u7jdtWg4CDtkvA/K7ADncwPgZcUrFOvAT4O/d8EvD3gPdlufV2CojtmEqxVBVvAdAgYNrmin+OKvbXbODugNe98H5dJ1exbI7bXsUvwKnA4wHzTwG+ds9HAnMC5glezSowwawHJGCZeRX/0FVs+0XgGvf8L3g1rW6VlunmPutxVPGrudKyt3LgCWYo8H2l9d0ETAlY94cHsP3vA14LXs00sOZwOLDGPZ9c6XvrQegJ5hzgf5Xi+Tfw5xC+1/OAz6r5TOcAH7vnyXg13yFBPn/gj6dWeIV7YA3tPOB993w6MAFoX8W6FDgp4PXlwLvu+bvA5QHzevLT/9X2AfPnAee65x/i1TabV9reDbgkFTDtv8CoYN95rDzsHEzdna6qOXi/dK4EPhCR1mFcf0fgBRHZKSI78RJOGd4/SYW1FU9UdRVeYdEf79f2q8B6EemJ9+v8A7doW7yqfMX78vB+GbWrar1BbNP92/bz8ZJVdQLX+R1eraC5iCSLyN0iskpEduMlb/BqJxU2VrOdtuy/D7SK2Ne56YHbbgsgIieLyBwR2e728SkB270H71fpWyKyWkRudNtYCfwWr/DaLCLPiEjbIJ/7QHUE2lZ8/y6+P1LN918Lge9pgVcDXxiwjTfddKi0fwk4bkLQERhaKf7z8Zq5KlT3vXbAq4FX5SWgl4h0xvtVv0tV5wWJIzD+jnjH3YaAmP6NV5MBuB4v6c4TkaUiMjbIuvYdR1T6n3LPG7D/d1XdZx2Hl7i/FpH5InJqQKxnVdp/R+L90Ix5lmDCRFXLVPV5vML/SDd5L94/boWaEo9WMW0tcLKq5gQ80lV1XZD3fYDXXJTqlvsAr+23CbDYLbMe7+AFQEQygWZ4TVjB4qmrwF52B+H9wtuK15w3Aq9GkI33iw+8f/SabAhcr4hIpe0AtHPTA7e9XkTS8Jrq/gG0cj8WXq/YrqruUdXrVLULcBrwOxE51s17SlWPxNuPCvwthFhrUnmfr8WrSQR+/41U9ZQg76ntdrbi1UR7B2wjW73OK1Bp/+Ltu0DBjvO1wAeV4s9S1ctCiHEtXnPTT4NXLcQ7H3kBcCEwo4Z1BX7etXg1mOYBMTVW1d5u3RtV9WJVbYvXAvBIpd6TlffFevd8v/8pN68Ur3kxeHCqK1T1PLwk9zfgWfc/uRavBhO4/zJV9e6a1hkLLMGEiXhG4BXiy9zkxcAZItLQHaDjKr1tE/v/A20CmolIdsC0x4A7RKSj204Lt51gPsCrTX3oXs92rz9S1TI37WlgjIj0d4XsncBcVf02yHorx3sgLhCRXiLSEK/56VkXUyO8f/pteIXVnbVY52tAbxE5w/UQupqfJvOWwNUikiIiZwGH4CWSVLxzPluAUhE5GdjXhVhEThWRbi457cL7AVEuIj1F5Bi37wr5sYNGXVXex/OAPSJyg7ueJVlE+ojI4DBsCwBVLcc7N3C/iLQEEJF2InKiW2QWMDrge/tzpVUEO85fBXqIyIVu36eIyGAROSSE0F4F2ojIb0UkTUQaicjQgPnT8ZrjTqPmBBP4eTcAbwH3ikhjEUkSka4icrT77GeJSHu3+A685BT43f5BRJq4SxKuAWa66U8D14pIZxHJwjuGZ2oIvQdF5AIRaeG+i51ucjlex4RfisiJ7rtPd5cztK9+bbHDEkzdvSIieXi9P+7Aaxtd6ubdj3eOYRPe+Y8nK733VmCaq/qerapf4x2kq920tsC/8E6+viUie/BO+A8luA/wCuyKBPMRXqFd8RpVfQf4P7xf7xvwzt2cW8N694u3hmWrMwOv3X0jXrPi1W76dLwmhXXAV3ifMySquhU4C7gbL0F1x+slFWium74V73s6U1W3qeoeF8MsvMLkN3j7u0J34B28c22fAo+o6vt4Selut76NeAnsplBjDmISXtPPThF50SXfU/GaPNe47T2OV8sLpxvwmgLnuCbKd/DOIaCqb+B1unjPLfNepfdWe5y7/XsC3rG1Hm9f/Q1v/wXl3ns88Ev3vhXA8ID5H+MVwotUtTbNduCdt0vFO9Z24PUyq2h2GgzMdf/XL+Odj1sd8N6X8DpeLMb7cTPJTZ+Md3x/iPddFeJ1pAjFScBSt81/4Z2bKVDVtXg1+z/i/Qhai9ejLy7Kbtm/WdqY+kdERuOdhD6ypmVNaEREge7uXJSfcbwHPKWqj0dpezHxueOFXXBkjIlLrplwIN4vfBOD4qKaZYwxgURkGl4z3m9dU5qJQdZEZowxJiKsBmOMMSYi6vU5mJycHO3WLfYH/927dy+ZmZl+h1EjizO8LM7wioc44yFGgIULF25V1RY1LxlcvU4wrVq1YsGCBX6HUaPZs2czbNgwv8OokcUZXhZneMVDnPEQI4CI1Lbbd5WsicwYY0xEWIIxxhgTEZZgjDHGRIQlGGOMMRFhCcYYY0xEWIIxxhgTEZZgjDHGRIQlGGOMMRFhCcYYY0xEWIIxxhgTEZZgjDHGRIQlGGOMMRFhCcYYY0xEWIIxxhgTEZZgjDHGRETcJRgROUlEvhGRlSJyo9/xGGOMqVpcJRgRSQYeBk4GegHniUgvf6MyxhhTlXi7o+UQYKWqrgYQkWeAEcBX4drAqi15vP7FBrbkFVFQXIaGa8VBbNxYxKtbPt/3WirNl4AJUmmuSPXLVl7bT5YNsp3K2xKBdeuKeG/Xkmre+5MNh7TemuL66eeTapetmLD2+2I+LVhWq+1Unp+SnER6SjIZ7pGWkkRGSjJZaQ1olpVG86xUmjRMJSmp+s9tTKIT1WgUoeEhImcCJ6nqRe71hcBQVb0yYJnxwHiAFi1aDJo1a1bI6/9+dxl/nVNIaTlkpkBKkhCN8qNcy0kSrzIZ7Nuo/FXV9M0Fzv/p16xBXu0/oeKpqiIiNW+3iveGsmzwCIPHuP9LRSunj1rsOwXKQ/i3ECA7TWidKbTOTKJtZhJdc5Lo2DiJBiEcOHl5eWRlZdW8IZ9ZnOETDzECDB8+fKGq5tZ1PfFWg6mRqk4AJgD07NlTa3P/6/HTF5CZXsZrVx9Fu5yMCEX4U/Fyn+5EirO8XCksLaOwpJyCkjIKissoLCljT2Ep2/YWsXVPEVvzitm4u5A1W/fy2eY8Zq8tBiA9JYnBnZpyYu/WnNi7NS0apUUszmiwOMMnHmIMp3hLMOuADgGv27tpdVZaVs4nq7Yxon/bqCYXE5uSkoSGqQ1omBra8qrKlj1FLPhuB/PWbOfD5Vv404tLuOWlJZzQqzWjf9aJoZ2bBm1KNKa+ibcEMx/oLiKd8RLLucBvwrHiZRv2kFdUytAuzcKxOpNgRISWjdM5pW8bTunbBlVl+aY8XvhsHc/M/543l26kf4ccbjr5YDvGTMKIq15kqloKXAn8F1gGzFLVpeFY96oteQD0atMoHKszCU5E6Nm6ETeefDBzbjqWO37Vh427CjlnwhwumbGAzbsL/Q7RmIiLtxoMqvo68Hq41/vDjnwA2jdpGO5VmwSXnpLM+UM78uuB7Zn00RoeeHcFx933Aed0T+Jo13HCmPoormowkbR2ewEtGqWRnpLsdyimnkpPSeaK4d1445qj6Nm6ERO/LOYPz35BYUmZ36EZExGWYJwfdubTvomd3DeR16VFFjPHH86Irik8t+gHfvXIJ6zdnu93WMaEnSUYZ9PuIlo3Tvc7DJMgkpKEX3VPZfKowazbkc8Zj37C0vW7/A7LmLCyBOPszC+mSWaIfVKNCZPhB7fkucuOICVJOOffc/h45Va/QzImbCzB4F3DsCO/hCYNU/wOxSSg7q0a8dzlR9AuJ4MxU+bz7rJNfodkTFhYggF2F5ZSVq40CfWqOmPCrE12BrMuOZyD2zTisicWWZIx9YIlGLzmMcASjPFVdsMUZowdui/JvPe1JRkT3yzBADvySwBokmlNZMZf2Q1TmDHOSzKXPrGIT1dt8zskYw6YJRhgx16vBpNjNRgTA7IzUpg+dggdmzbk4ukLWLLOepeZ+GQJBtjhmsiaWoIxMSKnYSrTxw0hOyOF0VPm8e3WvX6HZEytWYIBtu+1czAm9rTJzmD6uCGUK1w4ea6NX2bijiUYYGd+CUkCjdLjbmg2U891bZHFlNGD2ZZXzJip89lbVOp3SMaEzBIMsLOgmOyMFLv9rYlJ/Trk8PBvBrJsw26ueeYzykK53aYxMcASDLCnsJTGGdaDzMSu4Qe35LbTevPOss3c/tpXfodjTEisTQjIKywlK812hYltFx7eiTVb85n88Ro6Nctk1BGd/A7JmKCsVAX2FFmCMfHh5l8cwvfb87ntlaV0aJrBMQe38jskY6plTWR4NZhG6dZEZmJfcpLwwHn96dW2MVc+9ZmNwGximiUYYE9RifUgM3GjYWoDJo8aTE5GChdNW8DmPdZ92cQmSzDYORgTf1o2TmfCyFx25pdwyYyFdldME5MSPsGoKnlFpWRZDcbEmT7tsrn/nH589v1Obnr+S1St+7KJLQmfYIpKyykpU2siM3HppD5tuO74Hrzw2Toe/WCV3+EYs5+EL1Xz3JXRjayJzMSpK4/pxvLNedzz32/o1iKLE3q39jskYwCrwbCn0Esw1kRm4pWIcM+Zh3Jou2x+O3Mxyzbs9jskYwBLMORVJJg066Zs4ld6SjITRubSON3rWbY1r8jvkIyxBLOnyLvZmJ2DMfGuVeN0Jo7MZdveIi6dsZCiUutZZvyV8AnmxxqMJRgT//q2z+YfZ/VjwXc7uPmFJdazzPgq4UvVinMwVoMx9cWph7ZlxaY8/vXuCnq2asTFP+/id0gmQSV8qVrRi8xqMKY+uebY7qzYvIc731hG15aZNmaZ8YU1kRVZLzJT/yQlCfee1Z/ebRtz9dOLWb5pj98hmQSU8AlmT2EpqQ2SSGuQ7HcoxoRVRmoyE0fmkpGazLhp8/fdGtyYaEn4BJNXVGIXWZp6q012BhMuHMSm3UVc+sRCikvL/Q7JJJCYSzAico+IfC0iX4jICyKSEzDvJhFZKSLfiMiJ4djenkIbh8zUbwMOasI9Zx7KvDXbueUl61lmoifmEgzwNtBHVQ8FlgM3AYhIL+BcoDdwEvCIiNS5Xcu7F4wlGFO/jejfjiuGd+WZ+WuZ9NEav8MxCSLmEoyqvqWqpe7lHKC9ez4CeEZVi1R1DbASGFLX7dndLE2iuO74npzUuzV3vL6Md5dt8jsckwAklqvLIvIKMFNVnxCRh4A5qvqEmzcJeENVn630nvHAeIAWLVoMmjVrVtBt3PJxAc0yhGsGpkfkM4QiLy+PrKws37YfKoszvPyIs6hUuWteIRv2lnPz0HQOalxzI4Dtz/CJhxgBhg8fvlBVc+u6Hl9+uovIO0BVQ77erKovuWVuBkqBJ2uzblWdAEwA6Nmzpw4bNiz48vPeo1O7pgwb1r82mwmr2bNnU1OcscDiDC+/4uw/pJARD33MY0vhxSsPo2Wj4D+ubH+GTzzEGE6+NJGp6nGq2qeKR0VyGQ2cCpyvP1ax1gEdAlbT3k2rk/yiMjLTrIuySRytGqfz+KhcduSXcPF0uxumiZyYOwcjIicB1wOnqWp+wKyXgXNFJE1EOgPdgXl13d7e4lIyU+0cjEksfdpl889z+/PFDzu57j+fU14eu03lJn7FXIIBHgIaAW+LyGIReQxAVZcCs4CvgDeBK1S1Tj+9ysqVwpJyGlqCMQnoxN6tufGkg3ntiw38853lfodj6qGYK1lVtVuQeXcAd4RrW/nFXmc1ayIziWr8z7uwakseD7y3ki4tsjh9QDu/QzL1SCzWYKImv9irAFkNxiQqEeH20/sytHNTrn/2CxZ8u93vkEw9ktAJZm+R1WCMSW2QxGMXDKJtTjqXzFjI2u35Nb/JmBAkdIKxGowxniaZqUwaPZiSsnLGTp3P7sISv0My9UBCJ5h9NZhUq8EY07VFFo9dMIg1W/dy5VOfUVpmA2OauknoBLOvBmNDxRgDwBHdmnP76X34cPkW/vrqV36HY+JcQpese4utBmNMZecOOYhVW/KY+L81dG6eSSe/AzJxK7FrMEVWgzGmKjeefAjHHdKKv7z6FZ9tLq35DcZUIaETjNVgjKlacpLwwHn96dMum0c/L+KLH3b6HZKJQ9X+dBeRl0N4/3ZVHR2+cKLLepEZU72GqQ14fFQuJ9/7HmOnLuDFK46gfZOGfodl4kiwkvUQ4KIg8wV4OLzhRNfeolJSkoXUBgldkTOmWi0bpXNdbjp3LyhhzJT5PHvZEWRnpPgdlokTwRLMzar6QbA3i8htYY4nqvKLy6z2YkwN2mYl8diFgxg1eR6XzljItLFD7EeZCUm1R4mqBr9TV4jLxLK9RaV2/sWYEBzRtTl/+/WhfLp6Gzc+/wWxfKNCEzuqTTAi0lxE/iwiV4tIlog8KiJLROQlEal2QMp4kl9cZj3IjAnRGQPb87vje/D8onX8850Vfodj4kCweu5TQBo/3ndlNXAm8CrweORDizzvXjBWgzEmVFcd040zB7XnX++u4NmFP/gdjolxwX6+t1LVP4qIAN+p6j1u+tcickUUYou4/CI7B2NMbYgId53Rl427CrnxuS9ok53Oz7o19zssE6OC1WDKANwti7dWmlcvBinaW1xqIykbU0spyUk8csFAurTI5NInFrJ80x6/QzIxKliC6SIiL4vIKwHPK153jlJ8EWW9yIw5MI3TU5gyZggZKcmMmTKfzbsL/Q7JxKBgpeuIgOf/qDSv8uu4tLfIajDGHKh2ORlMHj2Ys//9KWOnzWfm+MPJtE4zJkC1R0NN18DUB1aDMaZu+rTL5qHfDOCiaQu44qlFTByZS0qyXSNjPMG6KX8pIl9U94hmkJGgqtaLzJgwOObgVtzxq77M/mYLN7/wpV0jY/YJ9vP9VPe3osfYDPf3AiDuj6DCknJUbSRlY8LhvCEHsXFXIf96dwWtszP43fE9/A7JxIBgTWTfAYjI8ao6IGDWDSKyCLgx0sFFko2kbEx4/fa47mzcVcgD766gVeM0zh/a0e+QjM9CaSwVEflZwIsjQnxfTNt3Lxg7B2NMWIgId/yqD8N7tuD/XlzC219t8jsk47NQEsU44BER+VZEvgUeAcZGNKoo2FeDsV5kxoRNg+QkHj5/IH3bZXPV04tY+N0Ov0MyPqoxwajqQlXtB/QD+qlqf1VdFPnQIivfJRirwRgTXg1TGzBp9GBaN07nomnzWbUlz++QjE+C9SI7NfC1qu5S1V3Blokne10TmdVgjAm/5llpTBs7hCQRRk2ex+Y9diFmIgpWg7lHRAaIyMDqHsCd0Qo03KwGY0xkdWyWyZQxg9m+t5gxU+aTV1Tqd0gmyoKVrpuA+2p4f9yO2b2vBmMJxpiIObR9Dg+fP5CLpi3gsicWMmnUYLtZWQIJ1k15WBTjiLp9NRhrIjMmoob3bMldZ/Tl+me/4MbnvuDes/vhDdJu6ruE/fm+t7iim7IlGGMi7ezcDmzaVci9by+ndXY61590sN8hmSiI2bqqiFwnIioizd1rEZEHRGSlG65mYF3WX+ASTHoDSzDGRMOVx3TjN0MP4pHZq5j2ybd+h2OiICZrMCLSATgB+D5g8sl4d9fsDgwFHnV/D0hBSRnpKUkkJVlV3ZhoEBH+clpvNu8u4tZXltI8K41fHNrG77BMBNVYgxGRhiLyfyIy0b3uHoXuyfcD17P/mGcjgOnqmQPkiMgBH50FxWVkpFjtxZhoapCcxIPnDWDQQU24duZiPllZ+V6Gpj6RmkY+FZGZwEJgpKr2EZGGwCeq2j8iAYmMAI5R1WvcyAG5qrpVRF4F7lbVj9xy7wI3qOqCSu8fD4wHaNGixaBZs2ZVuZ3Hvyziq21l3DesYSQ+Rq3k5eWRlZXldxg1sjjDK5Hj3Fui3Dm3gG0Fyo1D0umUXfcfe/GwP+MhRoDhw4cvVNXcOq9IVYM+gAXu72cB0z6v6X01rPMdYEkVjxHAXCDbLfct0Nw9fxU4MmAd7+Iln2q306NHD63O5U8u1OH/eL/a+dH0/vvv+x1CSCzO8Er0ODfsLNAj7npXB/7lLV29Ja/O64uH/RkPMaqqVpT7dX2EcpK/WEQycM1VItIVKKpjUjtOVftUfgCr8W7H/LmrvbQHFolIa2Ad0CFgNe3dtANSWFxmPciM8VHr7HSmjxuCAhdOmmu3Xa6HQkkwfwbeBDqIyJN4NYfrIxGMqn6pqi1VtZOqdgJ+AAaq6kbgZWCk6012GLBLVTcc6LYKSuwcjDF+69oiiymjvav9R06ex66CEr9DMmEUymCXbwNnAKOBp/GapWZHNqwqvY5Xw1kJTAQur8vK8ovLSLcEY4zv+nXI4bELBrFqSx4XT19AYUmZ3yGZMAk22GXgmGMdgQ3AeuCgul6DEipXk9nqnquqXqGqXVW1r1Y6uV9bhVaDMSZm/LxHC/5xVj/mrdnO1U9/RmlZud8hmTAIdh3Mve5vOpALfA4IcCiwADg8sqFFVkFJGRl2DsaYmDGifzu27y3mtle+4uYXlnD3r/vakDJxLthYZMMBROR5vPMgX7rXfYBboxJdBBXYSX5jYs6Yn3VmW14xD72/kuaNUvnDiTakTDwL5Ur+nhXJBUBVl4jIIRGMKSq8K/ktwRgTa647oQfb9hbx8PuraJ6VxpifdfY7JHOAQkkwX4jI48AT7vX5wBeRCyk67ByMMbFJRPjriD77msuaZqYyon87v8MyByCUbspjgKXANe7xlZsWt0rKyikpU0swxsSoBslJ/OvcAQzt3JTrZn3OB8u3+B2SOQChdFMuVNX7VfVX7nG/qsb1FVEFrhukneQ3JnalpyQzcVQu3Vs14rInFrJ47U6/QzK1FMpgl2tEZHXlRzSCi5TCYkswxsSDxukpTBszmGZZqYyZMo+Vm/P8DsnUQihNZLnAYPc4CniAH8/HxKV9NRhrIjMm5rVsnM6MsUNJThJGTprL+p0FfodkQhRKE9m2gMc6Vf0n8IsoxBYxlmCMiS+dmmcydcwQ9hSWcuGkuWzfW+x3SCYEoTSRDQx45IrIpcTojcpCte9ultZEZkzc6NMum8dH5fLDjgJGT5nHnkIbtyzWhdJEdm/A4y5gIHB2JIOKtIoEYzUYY+LL0C7NeOT8gSxdv5vx0xfauGUxLpQEM05Vh7vH8ao6Hojr+mlFE5ldyW9M/Dn2kFbce1Y/Pl29jats3LKYFkqCeTbEaXHDzsEYE99OH9CO207rzdtfbeKG576kvDz4nXmNP6o9lyIiBwO9gWwROSNgVmO8ATDj1r5zMJZgjIlbo47oxM78Eu5/Zzk5DVM4MtOSTKwJdrK+J3AqkAP8MmD6HuDiSAYVaYV2oaUx9cLVx3ZjR34xkz5aw47uKQwf7ndEJlCw0ZRfAl4SkcNV9dMoxhRx1kRmTP0gItxyai92F5Tw/GfrGDDnOy48rKPfYRknWBPZ9ar6dwIg/3QAABkgSURBVOA3InJe5fmqenVEI4ugfOtFZky9kZQk/O3MQ1n9w0ZueWkJjdMb2OCYMSJYE9ky97dOd46MRQUlZaQ1SCIpyW5mZEx9kJKcxOX905i0Mp3rZn1O4/QUhh/c0u+wEl6wJrJX3N9p0QsnOgqL7W6WxtQ3qcnC46NyOW/iHC57ciEzxg1lcKemfoeV0II1kb0CVNstQ1VPi0hEUVBg94Ixpl5qlJ7CtDFDOOvfnzJ26nxmjj+cXm0b+x1WwgrWRPaPqEURZQUl5ZZgjKmnmmWlMWPcUM589BNGTp7Hfy49nM7NM/0OKyFVe6Glqn5Q8QA+BXYA24FP3bS4VWBNZMbUa+1yMpgxbijlqlzw+Fw27orrW1jFrVAGu/wFsApvmP6HgJUicnKkA4ukgpJSq8EYU891a5nFtDFD2FVQwoWT5rLDRmCOulAHuxyuqsNU9WhgOHB/ZMOKrILiMruK35gE0Ld9NhNH5vLd9nxGT53P3qJSv0NKKKEkmD2qujLg9Wq8q/njVlFpOekpoXx0Y0y8O7xrMx7+zUCWrNvF+BkLKCq1EZijJZRSdoGIvC4io0VkFPAKMF9Ezqg0RlncKCotJ62B1WCMSRTH92rF3399KB+v3MY1Ty+2EZijJJQEkw5sAo4GhgFbgAy88clOjVhkEVRYUkaa1WCMSSi/HtSeW07txZtLN/LHF75E1QbHjLQa70ypqmOiEUg0WQ3GmMQ09sjO7Cwo4YF3V5DTMJWbTj4YERvRI1JqTDAi0hm4CugUuHw8X2hZWFJm52CMSVDXHtedXfnFTPhwNTkNU7h8WDe/Q6q3akwwwIvAJLxzL/Wi4dJqMMYkLhHhz7/sza6CEv7+5jdkZ6Rw/lAbgTkSQkkwhar6QMQjiZLycqXYepEZk9CSkoR7zurHnsJS/vTiEhqlp3Bav7Z+h1XvhFLK/ktE/iwih4vIwIpHJIMSkatE5GsRWSoifw+YfpOIrBSRb0TkxANZd7HrPWI1GGMSW0pyEg+fP5DBnZryu5mLef/rzX6HVO+EUoPpC1wIHMOPTWTqXoediAwHRgD9VLVIRFq66b2Ac/Fu49wWeEdEeqhqrTq1V9zN0mowxpj0lGQeH5XLb2wE5ogIpZQ9C+iiqker6nD3iEhycS4D7lbVIgBVrfhZMQJ4RlWLVHUNsBIYUtuVF5VaDcYY86PGbgTmtjkZjJ06n6Xrd/kdUr0hNfUFF5EXgfEBBX1kAxJZDLwEnAQUAr9X1fki8hAwR1WfcMtNAt5Q1WcrvX88MB6gRYsWg2bNmrXf+jfnl3P9hwVc3DeVn7VLifwHCkFeXh5ZWVl+h1EjizO8LM7wqmuc2wrKuWNuIaXlyh+HZtA6M/ytHPGyL4cPH75QVXPrup5QmshygK9FZD5QVDGxLt2UReQdoHUVs252MTUFDgMGA7NEpEuo61bVCcAEgJ49e+qwYcP2m7980x748EP69+3DsEPbHNgHCLPZs2dTOc5YZHGGl8UZXuGIs39uHmc99ikPfgnPXjaUNtkZ4QnOiZd9GS6hJJg/h3ujqnpcdfNE5DLgefWqVvNEpBxoDqwDOgQs2t5Nq5WKczBpDewcjDFmf11bZDF97BDOnTCHCyfNY9Ylh9M0M9XvsOJWjaVs4H1h3H1gyoCzIxjTi3gjNiMiPYBUYCvwMnCuiKS5iz+7A/Nqu/KKczA2mrIxpip92mXz+Khc1m7PZ/SUeeTZCMwHLKSf8SIyQETuEZFvgb8CyyIY02Sgi4gsAZ4BRqlnKTAL+Ap4E7iitj3IAIpK3El+60VmjKnGYV2a8cj5A1m6fjcXT1uwr+XD1E61payI9HDXv3wNPAh8j9cpYLiqPhSpgFS1WFUvUNU+qjpQVd8LmHeHqnZV1Z6q+saBrH9fN2XrRWaMCeLYQ1px71n9+HT1Nq56+jMbgfkABPsZ/zXetS6nquqRqvogXvNYXNvXTdlqMMaYGpw+oB23ndabt7/axA3PfUl5uY3AXBvBTvKfgXdh4/si8iZec1XcDztqNRhjTG2MOqITuwpKuO/t5TTOaMAtp/ayEZhDVG2CUdUXgRdFJBPvIsffAi1F5FHgBVV9K0oxhpXVYIwxtXXVMd3YmV/C5I/X0KRhKlcf293vkOJCKL3I9qrqU6r6S7yuwZ8BN0Q8sgixbsrGmNoSEf70i0P49cD23Pf2cqZ98q3fIcWFUK6D2UdVd+BdxDghMuFEnnVTNsYciKQk4W+/7svuwhL+/PJSsjNSOH1AO7/DimkJ9zO+ogaTmpxwH90YU0cNkpN48LwBHN6lGdf953Pe+WqT3yHFtIQrZYtKy0ltkERSkp2kM8bUXnpKMhNH5dKnbWOueGoRc1Zv8zukmJVwCaawpMzOvxhj6iQrrQFTxgyhQ9OGXDRtAV/+YCMwVyXhStqi0nI7/2KMqbOmmanMGDeE7IwURk2Zx8rNeX6HFHMSL8FYDcYYEyZtsjN44qKhJAlcOGku63YW+B1STEm4ktZqMMaYcOrcPJPpY4eSV1TKhY/PZWteUc1vShAJl2DsHIwxJtx6tW3MlNGDWb+rgFGT57G7sMTvkGJCwpW0RaXllmCMMWGX26kpj10wiOWb9nDRVBuBGRIywZRZE5kxJiKG9WzJfWf3Z/5327n8yUWUJPgIzAmXYApLrAZjjImcX/Zry+2n9+G9rzfz+/98ntAjMNdqqJj6wGowxphIO39oR3YVlPD3N78hOyOF207rnZAjMCdcgrEajDEmGi47uiu78kv494eryclI4Xcn9PQ7pKhLuARjNRhjTDSICDeefDA780t44L2VNM5IoZvfQUVZwv2UtxqMMSZaRIQ7z+jLKX1bc/try/jfD4nVfTnhStqi0jLSrAZjjImS5CTh/nP6c1T35kxeUsx/l270O6SoSagEo6oUlpSTbjUYY0wUpTVI5rELBtElO4mrnvqMj1du9TukqEiokra4rOJ2yVaDMcZEV2ZaA64dlE7n5plcPH0Bi9fu9DukiEuoBFNY4hKM1WCMMT7IShVmjBtC86w0Rk+Zx/JNe/wOKaISqqQtLrUEY4zxV8vG6TwxbiipyUlcOGkua7fn+x1SxCRUSVvRRJZqCcYY46ODmjVkxrihFJaUc8GkuWzeU+h3SBGRUCVtRQ3GEowxxm89WzdiypjBbN5dxMhJ89iVX/+6MCdUSVsx8FxKckJ9bGNMjBp4UBMmjBzEqi15jJ02n/ziUr9DCquEKmn31WAswRhjYsRR3Vvwr3MH8Nn3O7j0iUX7yqn6IKFK2iJrIjPGxKBT+rbhrjP68uHyLVw7azFl9WQE5oQai6yiicxqMMaYWHPO4IPYVVDCna9/TeP0FO78VZ+4H4E55kpaEekvInNEZLGILBCRIW66iMgDIrJSRL4QkYG1Xbed5DfGxLLxP+/K5cO68vS87/n7f7/xO5w6i8UazN+B21T1DRE5xb0eBpwMdHePocCj7m/ILMEYY2LdH07syc6CEh6dvYqcjBQuObqr3yEdsFhMMAo0ds+zgfXu+QhguqoqMEdEckSkjapuCHXFJXYdjDEmxokIfx3Rh90FJdz1xtc0zUzlrNwOfod1QMQrr2OHiBwC/BcQvCa8I1T1OxF5FbhbVT9yy70L3KCqCyq9fzwwHqBFixaDZs2atW/ep+tL+fcXRdx1ZAZtsmInyeTl5ZGVleV3GDWyOMPL4gyveIizNjGWlCv/XFjIsu3lXD0gjf4to1cfGD58+EJVza3zilQ16g/gHWBJFY8RwAPAr91yZwPvuOevAkcGrONdIDfYdnr06KGBZs7/Xjve8Kp+v22vxpL333/f7xBCYnGGl8UZXvEQZ21j3FNYoqc+8D/t+afXdcG32yMTVBWABRqGst6Xn/Gqepyq9qni8RIwCnjeLfofYIh7vg4IrCe2d9NCZmORGWPiSVZaA6aMGUzrxumMnTqfFXE2OGYslrTrgaPd82OAFe75y8BI15vsMGCX1uL8C9iV/MaY+NM8K43pY4eSkpzEyMnzWL+zwO+QQhaLJe3FwL0i8jlwJ+58CvA6sBpYCUwELq/tiq0XmTEmHh3UrCHTxg5mT2EpoybPY2d+sd8hhSTmSlpV/UhVB6lqP1UdqqoL3XRV1StUtauq9tVKJ/dDYQnGGBOverfNZsLIQXy3LZ9x0xZQUFzmd0g1SqiStqKJrEFSfF8da4xJTEd0bc4/z+3Pou93cNXTiygti+1xyxIqwRSVlZPaICnuh18wxiSuU/q24S8j+vDOss388YUvK3rVxqRYvNAyYopLy0mzE/zGmDh34WEd2bKniAfeXUGLRmn84cSD/Q6pSgmXYOz8izGmPrj2uO5s2VPEw++vonlWGmN+1tnvkH4ioRJMSVm5dVE2xtQL3pAyvdmWV8RfXv2KZllpnNavrd9h7SehSlurwRhj6pMGyUk8cN4ABndsynWzFvPRiq1+h7SfhCpti8sswRhj6pf0lGQmjsqla4ssLpmxgC9/2OV3SPskVGlbXKrWRGaMqXeyM1KYNnYIOQ1TGT1lHmu27vU7JCDREozVYIwx9VSrxulMHzeEclVGTp7Llj1FfoeUYAmmtMy6KRtj6q2uLbKYMmYIW/YUMXbqfPYWlfoaT0KVtnaS3xhT3/XvkMND5w1k6fpdXPmUv1f7J1RpW1KmpCTbVfzGmPrtuF6t+OvpfXj/my386cUlvl3tn1DXwVgNxhiTKM4f2pH1Owt4+P1VtM3J4Opju0c9hsRKMGXlpDZI9jsMY4yJit+f0JMNuwq57+3ltMlO56zcDjW/KYwSK8GUllsTmTEmYYgId59xKJt3F3HT81/SsnE6R/doEbXtJ1R7UXFZud0u2RiTUFIbJPHoBQPp3qoRlz+xkCXronchZkKVtsWl5aRaN2VjTIJplJ7C1DGDyc5IYczU+azdnh+V7SZUaWsn+Y0xiapV43SmjR1CUUkZo6dE57bLCVXa2mjKxphE1r1VIyaOzGXt9gIunr6AwpLI3nY5YUrb8nKltFytBmOMSWhDuzTjvnP6Mf/bHVw7czHl5ZG7RiZhSttidzWrJRhjTKI79dC2/OkXh/DGko3c/tqyiG0nYbop70sw1kRmjDFcdFQX1u8sZPLHa2ibk85FR3UJ+zYSJ8GUWg3GGGMC/ekXh7BxdwG3v7aMtjkZnNK3TVjXnzCl7b4EYzUYY4wBIClJuO/s/gzq2IRrZy5m4Xc7wrv+sK4thpXYORhjjPmJ9JRkJo7MpU12OhdPX8B328J3s7KEKW0rajDWTdkYY/bXNDOVKWOGoKqMnjI/bOtNmNK2yM7BGGNMtTo3z2TiyFzW7SwI2zoTprS1bsrGGBNcbqem3HtWv7CtL2FK2xI7yW+MMTX6Zb+2YVtXwpS2VoMxxpjoSpjS1ropG2NMdPlS2orIWSKyVETKRSS30rybRGSliHwjIicGTD/JTVspIjfWdpvWTdkYY6LLr9J2CXAG8GHgRBHpBZwL9AZOAh4RkWQRSQYeBk4GegHnuWVDVmTdlI0xJqp8GSpGVZeBdzvPSkYAz6hqEbBGRFYCQ9y8laq62r3vGbfsV8G2893ucnrf8iYAJW7EULujpTHGREesjUXWDpgT8PoHNw1gbaXpQ6tagYiMB8a7l0Vf/fXkJYHzD/pbeAINs+bAVr+DCIHFGV4WZ3jFQ5zxECNAz3CsJGIJRkTeAVpXMetmVX0pUttV1QnABBfDAlXNreEtvrM4w8viDC+LM3ziIUbw4gzHeiKWYFT1uAN42zqgQ8Dr9m4aQaYbY4yJQbF2QuJl4FwRSRORzkB3YB4wH+guIp1FJBWvI8DLPsZpjDGmBr6cgxGRXwEPAi2A10RksaqeqKpLRWQW3sn7UuAKVS1z77kS+C+QDExW1aUhbGpCZD5B2Fmc4WVxhpfFGT7xECOEKU5Rjdz9mI0xxiSuWGsiM8YYU09YgjHGGBMR9SLB1DSMjOs0MNPNnysinXyIsYOIvC8iX7lhcq6pYplhIrJLRBa7xy3RjtPF8a2IfOli+El3RfE84PbnFyIy0IcYewbsp8UisltEfltpGV/2p4hMFpHNIrIkYFpTEXlbRFa4v02qee8ot8wKERnlQ5z3iMjX7nt9QURyqnlv0GMkCnHeKiLrAr7bU6p5b52GmKpjjDMD4vtWRBZX895o7ssqy6GIHZ+qGtcPvJP+q4AuQCrwOdCr0jKXA4+55+cCM32Isw0w0D1vBCyvIs5hwKsxsE+/BZoHmX8K8AYgwGHA3Bg4BjYCHWNhfwI/BwYCSwKm/R240T2/EfhbFe9rCqx2f5u4502iHOcJQAP3/G9VxRnKMRKFOG8Ffh/CcRG0bIhkjJXm3wvcEgP7sspyKFLHZ32owQzBDSOjqsVAxTAygUYA09zzZ4FjpYpxaiJJVTeo6iL3fA+wjB9HKYg3I4Dp6pkD5IhIGx/jORZYparf+RjDPqr6IbC90uTAY3AacHoVbz0ReFtVt6vqDuBtvDH5ohanqr6lqqXu5Ry8a858Vc3+DEUoZUNYBIvRlTVnA09HYtu1EaQcisjxWR8STDt+OoxM5YJ73zLun2cX0Cwq0VXBNdENAOZWMftwEflcRN4Qkd5RDexHCrwlIgvFG3qnslD2eTSdS/X/vLGwPwFaqeoG93wj0KqKZWJtv47Fq6lWpaZjJBqudE15k6tp0omV/XkUsElVV1Qz35d9WakcisjxWR8STFwRkSzgOeC3qrq70uxFeM08/fCuE3ox2vE5R6rqQLzRq68QkZ/7FEeNxLvw9jTgP1XMjpX9uR/12hti+voAEbkZ71q0J6tZxO9j5FGgK9Af2IDXBBWrziN47SXq+zJYORTO47M+JJhgw8v8ZBkRaQBkA9uiEl0AEUnB+1KfVNXnK89X1d2qmueevw6kiEjzKIeJqq5zfzcDL/DjiNYVQtnn0XIysEhVN1WeESv709lU0Yzo/m6uYpmY2K8iMho4FTjfFTY/EcIxElGquklVy1S1HJhYzfZ935+uvDkDmFndMtHel9WUQxE5PutDggllGJmXgYoeD2cC71X3jxMprh12ErBMVe+rZpnWFeeGRGQI3vcT1UQoIpki0qjiOd5J3yWVFnsZGCmew4BdAdXraKv212Es7M8AgcfgKKCqAV//C5wgIk1ck88JblrUiMhJwPXAaaqaX80yoRwjEVXpnN+vqtl+LAwxdRzwtar+UNXMaO/LIOVQZI7PaPRciPQDr1fTcrweIze7aX/B+ycBSMdrQlmJN7ZZFx9iPBKv2vkFsNg9TgEuBS51y1wJLMXr7TIHOMKHOLu47X/uYqnYn4FxCt4N4FYBXwK5Pn3vmXgJIztgmu/7Ey/hbQBK8Nqpx+Gd83sXWAG8AzR1y+YCjwe8d6w7TlcCY3yIcyVeO3vFMVrR+7It8HqwYyTKcc5wx94XeIVjm8pxutc/KRuiFaObPrXieAxY1s99WV05FJHj04aKMcYYExH1oYnMGGNMDLIEY4wxJiIswRhjjIkISzDGGGMiwhKMMcaYiLAEYxKCiJTJ/qMvd/I7pnAQkdEiskVEHnevh4nIq5WWmSoiZwZZxz0islFEfh/peE1i8eWWycb4oEBV+1c1w118JupdFR6PZqrqlQf6ZlX9g4jsDWdAxoDVYEyCEpFO7j4h0/GunO4gIn8QkfluAMXbApa9WUSWi8hHIvJ0xS99EZktIrnueXMR+dY9T3a1gop1XeKmD3PveVa8e648GTDSwGAR+cQNzDlPRBqJyIci0j8gjo9EpF8dPnNuQA3uSxGxi+BMRFkNxiSKDPnxhk9rgGuB7sAoVZ0jIie410PwRip42Q06uBdviJH+eP8vi4CFNWxrHN7wOYNFJA34WETecvMGAL2B9cDHwM9EZB7eWFXnqOp8EWkMFOAN6TEa+K2I9ADSVfXzED7rUbL/za0OwrsvzgL3ORCRe4A3Q1iXMQfMEoxJFPs1kblzMN+pdz8b8MZVOgH4zL3Owks4jYAX1I3LJSKhjGV1AnBowHmPbLeuYmCeunGpXBLohHf7iA2qOh+8QTrd/P8A/ycif8AbomNqiJ/1f6p6asBn3e99InIO3s2xTghxfcYcEEswJpEFnncQ4C5V/XfgAlLpNsyVlPJjM3N6pXVdpar7DQQoIsOAooBJZQT5H1TVfBF5G+9mUGcDg4LEEhIR6YN3N8ifq2pZXddnTDB2DsYYz3+BseLdJwMRaSciLYEPgdNFJMONevvLgPd8y4+F/pmV1nWZGxYdEenhRsqtzjdAGxEZ7JZvJN4w7wCPAw8A89W7i+ABE5EcvEEZR6rqlrqsy5hQWA3GGLxbBYvIIcCn7rx7HnCBqi4SkZl4o91uxhsCvsI/gFni3YXwtYDpj+M1fS1yJ/G3UPUtaCu2XeyarR4UkQy88y/HAXmqulBEdgNTwvAxRwAdgYnuM1JdzzpjwsFGUzamFkTkVryC/x9R2l5bYDZwcFXdqMW7OVhuXbopu/XcShQ/l0kM1kRmTIwSkZF490u/Ocg1OgXAyRUXWh7gdu4BLmD/c1LG1JnVYIwxxkSE1WCMMcZEhCUYY4wxEWEJxhhjTERYgjHGGBMRlmCMMcZExP8DiPV7dOUyyIkAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "from matplotlib import pyplot as plt\n",
        "import numpy as np\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "ax.plot(w, 20 * np.log10(np.maximum(abs(h), 1e-5)))\n",
        "\n",
        "ax.set_title('Butterworth bandpass filter frequency response')\n",
        "ax.set_xlabel('Frequency [Hz]')\n",
        "ax.set_ylabel('Amplitude [dB]')\n",
        "ax.axis((0, 20, -100, 10))\n",
        "ax.grid(which='both',\n",
        "        axis='both')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G26eXYXnpHWF"
      },
      "source": [
        "<div class=\"alert alert-block alert-info\"><p><b>Question:</b> What does this plot tell us about the filter characteristics? What types of noise does the filter attenuate?</p></div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JJD2n5hGpHWG"
      },
      "source": [
        "<div class=\"alert alert-block alert-warning\"><p><b>Explanation:</b> This function generates the co-efficients for a Butterworth filter. The filter-type is specified as 'bp' - a bandpass filtter. The filter frequencies are specified in Hz (because the sampling frequency, fs, has also been specified): a high-pass frequency of 0.7 Hz, and a low-pass frequency of 10 Hz.</p></div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dnSasCJ5pHWG"
      },
      "source": [
        "<div class=\"alert alert-block alert-info\"><p><b>Extension 1:</b> How could we re-design the filter to retain frequency content of up to 20 Hz, but eliminate mains frequencies?</p></div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3cLPjh14pHWG"
      },
      "source": [
        "<div class=\"alert alert-block alert-info\"><p><b>Extension 2:</b> What would be appropriate cut-off frequencies when using the PPG for different purposes, e.g. heart rate monitoring, or blood pressure estimation? See <a href=\"http://peterhcharlton.github.io/publication/wearable_ppg_chapter/\">this book chapter</a> (Sections 2.2.4 to 2.2.5 on Sampling Frequency and Bandwidth) for details.</p></div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s8NgRIcapHWH"
      },
      "source": [
        "---\n",
        "## Filter the PPG signal"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5QTJpU4lpHWH"
      },
      "source": [
        "- Filter the PPG signal in preparation for differentiation."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "AHbccVHppHWH"
      },
      "outputs": [],
      "source": [
        "ppg_filt = sp.sosfiltfilt(sos_ppg, ppg)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kw90DmO_pHWH"
      },
      "source": [
        "- Plot original and filtered PPG signals"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "id": "g7E41qWdpHWI",
        "outputId": "40f0345d-83c9-439d-f449-72a30cce15c0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        }
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd3gUVffHvzeQkIA0IfQSRJqA0hUUBASkSBEVwQKIr4gKtteOBXtDf9gVfQULWEBAsIDSFFAwQZDeawDpPRBIcn5/nL3MJtkyM1tmkj2f59nnzs7cuXN3dnfOveece44iIgiCIAiCVeKc7oAgCIJQMBEBIgiCINhCBIggCIJgCxEggiAIgi1EgAiCIAi2KOp0B8JN+fLlKSUlxeluCIIgFCiWLl16gIiSrZxT6ARISkoK0tLSnO6GIAhCgUIptd3qOaLCEgRBEGwhAkQQBEGwhQgQQRAEwRaFzgbii7NnzyI9PR2nT592uisxSWJiIqpVq4b4+HinuyIIQhiJCQGSnp6OkiVLIiUlBUopp7sTUxARDh48iPT0dNSqVcvp7giCEEYcVWEppT5VSu1TSq3yc1wppd5WSm1SSq1QSjWzc53Tp0+jXLlyIjwcQCmFcuXKyexPEAohTttAxgPoGuB4NwB1PK+hAD6weyERHs4h914QCieOqrCI6HelVEqAKr0BfE4cc36xUqqMUqoyEe2JSgcFQRBcwtatwIgRwLFj/uvceitwxx3R65PTM5BgVAWw0+t9umdfLpRSQ5VSaUqptP3790etc5Gge/fuOHLkSMA6Tz/9NGbPnm2r/fnz5+Oaa67xub906dJo0qQJGjRogGeffTbgfgD466+/0L59e9SpUwfNmjVDjx49sHLlSlv9EgQhMOPHAz/+CCxY4P+1bVt0+1QojOhENBbAWABo0aJFgcyQRUQgIvz0009B6z733HMR6UPbtm3xww8/4OTJk2jSpAl69uzpd3/VqlXRr18/TJw4EW3atAEALFy4EJs3b0bjxo0j0j9BiGXWrOHymWeAjh1916lePXr9Adw/A9kFwPuWVPPsK3C8+eabaNSoERo1aoQxY8YAALZt24Z69eph4MCBaNSoEXbu3ImUlBQcOHAAAPD888+jXr16uOKKKzBgwACMHj0aADB48GBMnjwZAIdueeaZZ9CsWTM0btwY69atA8Czg9atW6Np06Zo06YN1q9fb7qvJUqUQPPmzbFp0ya/+999910MGjTonPAAgCuuuAJ9+vSxf5MEQfCLFiC9egHt2vl+RdvR0e0CZDqAgR5vrMsAHA3V/qFUZF6BWLp0KcaNG4clS5Zg8eLF+Pjjj7Fs2TIAwMaNG3H33Xdj9erVqFmz5rlzUlNT8d133+Gff/7Bzz//HDC+V/ny5fH333/jrrvuOidk6tevjwULFmDZsmV47rnn8MQTT5i+RwcPHsTixYvRsGFDv/tXr16NZs1sOcUJgmARImDjRt6uX9/ZvnjjqApLKfUVgPYAyiul0gE8AyAeAIjoQwA/AegOYBOADAC3OdPT0Fi4cCGuvfZalChRAgDQt29fLFiwAL169ULNmjVx2WWX5Ttn0aJF6N27NxITE5GYmHhOneSLvn37AgCaN2+OKVOmAACOHj2KQYMGYePGjVBK4ezZs0H7uWDBAjRt2hRxcXF47LHH0LBhQ8yfP9/n/rxceumlOHbsGLp06YK33nrL1H0RBMEcp04BZ88CSUlA8eJO98bAaS+sAUGOE4B7wnvNcLYWOlqohEKxYsUAAEWKFEFWVhYA4KmnnkKHDh0wdepUbNu2De3btw/ajrZ1mNnfsGFD/P333+jduzcAYMmSJZg8ebLP8wVBCA3teVWqlLP9yIvbVViFgrZt22LatGnIyMjAyZMnMXXqVLRt2zbgOZdffjlmzJiB06dP48SJE5YfzEePHkXVquywNn78eLtd98s999yD8ePH448//ji3LyMjI+zXEQTBvQKkUHhhuZ1mzZph8ODBaNWqFQDgP//5D5o2bYptAXzuWrZsiV69euHiiy9GxYoV0bhxY5QuXdr0NR955BEMGjQIL7zwAnr06BHqR8hHpUqV8M033+DRRx/Frl27UKFCBZQvXx5PP/102K8lCLGOWwWIIrfpdEKkRYsWlNfgvHbtWjRo0MChHtnnxIkTOO+885CRkYF27dph7NixBdZwXVC/A0FwA3PmAJ06AR06AHPnRuYaSqmlRNTCyjkyA3ExQ4cOxZo1a3D69GkMGjSowAoPQRBCw60zEBEgLmbixIlOd0EQBBfgVgEiRnRBEASXIwJEEARBsIUIEEEQBMEWIkAEQRAEW4gAiXHefvttNGjQADfffDOmT5+OV155BQAwatSoc/Grxo8fj927d0fk+tu2bUOjRo187k9KSkKTJk1w0UUXYdiwYcjJyfG7H+D4Xddccw1q166N5s2bo0OHDvj9998j0m9BENwrQMQLK0q8//77mD17NqpVqwYA6NWrV74648ePR6NGjVClShXT7WZlZaFo0dC+xtq1a2P58uXIyspCx44dMW3aNDRr1szn/u7du6NHjx4YPXr0uc+watUqpKWloV27diH1QxAE34gAiWGGDRuGLVu2oFu3bhgyZAjKli2LtLQ0vPvuu+fqTJ48GWlpabj55puRlJSEP//8E2vWrMGDDz6IEydOoHz58hg/fjwqV66M9u3bo0mTJli4cCEGDBiA9u3b+6y3dOlSDBkyBADQpUuXoP0sWrQo2rRpg02bNuVac+K9f8KECWjdunUuAajD1AuCEBncKkBiT4XlQDz3Dz/8EFWqVMG8efPwwAMP+Kxz/fXXo0WLFpgwYQKWL1+OokWLYsSIEZg8efI5QTBy5Mhz9c+cOYO0tDTce++9fuvddttteOedd/DPP/+YujUZGRmYM2dOvoRQ3vsljLsgRB+3ChCZgbiU9evXY9WqVejcuTMAIDs7G5UrVz53/MYbbwxY78iRIzhy5Mg5tdKtt96Kn3/+2ee1Nm/ejCZNmkAphd69e6Nbt27Ytm2bz/2//vprrnOvvfZabNy4EXXr1j0XSl4QhPAiAsQtFJDYX0SEhg0b4s8///R5XIeB91cvWF51b7Stw8z+hg0b5jKYT506FWlpaXjooYdMX08QBGu4VYDEngrLxZQsWRLHjx8HANSrVw/79+8/JxjOnj2L1atX5zvHX70yZcqgTJkyWLhwIQBgwoQJYenjTTfdhEWLFmH69Onn9kkYd0GILG4VILE3A3ExgwcPxrBhw84Z0SdPnox7770XR48eRVZWFu6///582QATEhL81hs3bhyGDBkCpZQpI7oZkpKS8MMPP+DBBx/E/fffj4oVK6JkyZJ48sknw9K+IAi5ycwEzpwB4uMBT+441yDh3IWoIN+BINhjzx6gShUgORnYty9y17ETzl1UWIIgCC7m0CEuy5Vzth++EAEiCILgYg4e5FIEiIMUNlVdQULuvSDYRwSIwyQmJuLgwYPyIHMAIsLBgweRmJjodFcEoUDiZgHiqBeWUqorgLcAFAHwCRG9kud4DQCfASjjqfMYEf1k9TrVqlVDeno69u/fH4ZeC1ZJTEw8FwNMEARriADxgVKqCID3AHQGkA4gVSk1nYjWeFV7EsC3RPSBUuoiAD8BSLF6rfj4eNSqVSsMvRYEQYgubhYgTqqwWgHYRERbiOgMgK8B9M5ThwDopTOlAUQm1rkgCIJLEQHim6oAdnq9T/fs82YUgFuUUung2ccIXw0ppYYqpdKUUmmiphIEoTAhAsQ+AwCMJ6JqALoD+EIpla/PRDSWiFoQUYvk5OSod1IQBCFSiADxzS4A1b3eV/Ps8+Z2AN8CABH9CSARQPmo9E4QBMEFiADxTSqAOkqpWkqpBAD9AUzPU2cHgKsAQCnVACxAREcVYxABy5cD27c73ROhIHPsGHDzzcDFFwOtWwMrVjjdI3O4WYA45oVFRFlKqeEAZoFddD8lotVKqecApBHRdAD/BfCxUuoBsEF9MMlijpjj/vuBt98GSpQA1q0DxCNY0Jw4AYwdCxw9CgwcCNSu7b/uf/8LTJxovL/0UmDbNqBixYh30zZE7g5lEhPBFIWCy7JlgHcCxKuv5vK664A77nCmT4I72L0b6NsXWLKE31euDKSmAlXzuuKAR/FVq3Jk23HjgP/8B8jOBm69Ffj88+j22wpHjwJlygDnnQd4Mj1EDAmmKBQ6vvySy65dgZIlgVmz+DV0qLP9Epzl5EngkktYeCQlARUqcNRaf7+Ln39m4dGpEzB4MPDJJ7z/66+Bv/+OWrct42b1FSACRHA5S5dyOWIE8PzzuY8dOBD9/gjuYOxY4/v/7Tdg5UogMRH46Sdg06b89X/5hcvu3bkcPJhnsGfPAo8/HpUu22LvXi7dqmYTASK4lpwcVmEBrMa64Ybcx0VTGZssWgQ8+ihvT5kCtGzJM5C+fXnfjBm56xMZAsQ7r9orr7DQ+eUXYFde/0+X8O+/XFaq5Gw//CECRHAt6ensOZOczH+gKlVyH09NdaZfgnP8+itwzTU8c7jvPuDaa41j3bpxOWtW7nNWruSRfJUqwEUXGfvPPx9o1463Fy+ObL/tIjMQQbDJxo1c1q1r7Js2DShdmrdFgMQW//sfzyCOHGFbxujRuY937szlb78Bp08b+xct4rJDB0Cp3Oe0asXlX39Fps+hIjMQQbCJFiB16hj7evc2BEdqKqsnhMLP6tXAXXfx9vXXA99+CxTNswihYkWgSRMWHgsWGPv170ULC29atsxdx21oASIzEEGwiC8BAgAXXsiujf/+617dtRBeRo9mtdWQIcCkSUDZsr7raTdvbfMAzAmQpUvZ5uY2tApLZiCCYBF/AkQpoIXHW92tI0chfMydC4wfD8TFAY88EriuNpJrAXLiBLBmDc9WmjTJX79yZV4fcuwYsGFDWLsdFkSFJQg28SdAAPerHoTwQAS8/DJvjxwJ1KsXuH6bNiwsVq7khXd//80zi8aN2ePKF3pm4sbfkhjRBcEG2dnAli28feGF+Y+LAIkNvvwSmD2bw9iM8JnMITeJiRzrioiFx2+/8f7LLvN/jv4tuc2QTiQzEEGwxY4dwJkz7Hp53nn5j+s/fVqaGNILM2+9xeWbb7I7txkuvZTLX34BZs7kbe/1H3nR6lC3rUg/fhw4dQooXtz3f8ANiAARXEkg9RXAeutKldil09fKY6Hgk5rKxu3zz+dAiWYZMIDLl14C/viDZyUdO/qvr9eGuM0G4nYDOiACRHApwQSIUqLGKuy8/z6Xt93m337hiyuuABo2NN737w+UKuW/fpUqrCI7cMCIfOsG3O7CC4gAEcLAH38AvXrx4i79uuuu0KKHBhMggAiQwsyJExzoEADuvNPauUoBDz5ovNfrRwLV14tV3TQLKQgzEMfygQiFg6VL2ff+xInc++fMAdauZT10QoL1drdu5TJQfgcRIIWX+fN5QWCrVoEHEf647TYO956VZfxOAlG3Lsdd27AhsME9mrjdgA7EuAAhYuOcXrV6ySXAU0/lX+Eq+GfkSBYe/fsDt9/O+06d4lHjb78BU6cCN95ovd30dC6rV/dfx9v4mZUl31thYronN2nXrvbOVwp48knz9bV78Pr19q4XCdzuwgvEuAB56aXcP7Lvv2eh8txzzvWpIHHkCM80ihQB3nkHKO+VrX7TJlYjTJoUmgDxlRxIU748UKsWz1bWrGH3TaHgc/KkkTnQzm/HDm5UYRWEGUjM2kB27QKefppHKmPGAM88w/vffZdXpQrB+eEHHvlfeWVu4QFwvKIiRTj4oV7PYZYzZ4B9+3jlcbDRl6ixCh+pqSxEmjXLHT03kmgBom1vbkCM6C5m5kxeodqjB4eFHjUKaN0aOHyYc3ALwfnuOy51HgZvqldntVZ2du481GbYvZvLypWDq6VEgBQ+9Hep13NEgxo1uNQzXzeg47zlTWPgJmJWgMydy6UOvgYAn37Ko94vvzT0j4JvTp40Fmn16eO7jvbHnzLFWtv6jxNIfaURAVL4+PVXLn0FP4wUyck8Yz54kFPfuoEdO7jUws2NxKwA0Q+cNm2MffXrG8lqdM5kwTfz5rGXzKWX+n/QX3UV5zFftszwqjKDFiDVqgWv26wZqyFXrAjNbVhwBzt3sgBJTGTX8GgRF8czXoBzqztNRgYLs/h4zrboVmJSgBw5wrrOYsU4yJo399zD5auvcmazpUvZQCvhMnLz++9cXnWV/zqJiawiBHh2ZxYzBnRNyZJA27ZsizGrKtuzh7/XpUt5JiW4h6VLubzySl6BHk30780NKQJ27uSyenUWbm7F0a4ppboqpdYrpTYppR7zU6efUmqNUmq1UsqiNt03q1Zx2agRS3hvOndmz4/jx9mFsEULXtX6zjvhuHLhYc4cLnVKUH8MH87l+++bz7dgZQYCAIMGcalVaoHYvp3XlrRowa86dYBt28xdR4g8a9Zw6b2SPFpoW4PbBIibcUyAKKWKAHgPQDcAFwEYoJS6KE+dOgAeB3A5ETUEEBbztvb1rl/fV784deZtt7F6pFYt3u/WnMlOsH07r70oUSK4ALn8ctbhHjpkCO5gWJmBAIYaMpgdZNs29rY5dYrVFdWr82ykTx+ewQjOowVItLyvvElJ4dKKujVSFBQB4uQ6kFYANhHRFgBQSn0NoDeANV517gDwHhEdBgAi2heOC2tfb3+5BUqUMFQuCxbwQ9KqK2ph5r33uLzmGiApKXj9tm2BCRN4YaGZtRr6D2zWeFi3Lsc62rWLv6cLLsh9/NAh4LXXgA8/ZBdhAPjxRx4cNGkC/PMPrwG67jpz13MTqanAxx+z3c4XlSuzh6GdaABOsG4dl74Gd5FGRz1wQ3DOgmBAB5wVIFUB7PR6nw4gr+NeXQBQSi0CUATAKCLKp6hQSg0FMBQAapi443oGEiw5DWDMQNwwKnEDBw8aQe4eesjcOV26sAD5/vvgOR2IrD9E4uKAnj35Gh9/bCQgAritG24wZj8XXsiCTKsr7r8feOAB4N57uZ8lS5q7ptMQsbfgsGFscA3Ev//yrFop49wlS3j21amTez4zkbkYaJFC553ZvDn6185LQZmBgIgceQG4HsAnXu9vBfBunjo/AJgKIB5ALbDAKROo3ebNm1MwGjQgAoiWLQtalbKziRISuP7y5cHrF3aefJLvRdeu5s85dIioaFGiIkWIDhwIXHfnTm6/XDlr/Vq0iM8rX57o9Gned+wYUdmyvL9sWaJ33yU6fDj3eZmZRM2acZ3337d2TSeZOpX7DBD17En06af5Xw88YNSZONE49+WXjf21a/M9z8lx7rNo9u7lPpUq5Ux/Nm7k69eoEf1r56VzZ+7Ljz9G75oA0sjqc9zqCeF6AWgNYJbX+8cBPJ6nzocAbvN6PwdAy0DtBhMgZ88SxcfzJz9xwtyN7d+f65cpQ/Tnn+bOKYykphrCdNEia+d26cLnffpp4HqzZ3O9K66w1n5ODlGTJnzuF1/wvpdeMh6Uq1b5P/err7hO48bueJCaoUMH7nPz5kQnT/quc+aMcU8Aot9+I5o2zXjv/erWzfnPrgcBJsaAEeHUKb5+0aI8cHSS+vW5LytWRO+aBU2AFAWwxTOzSADwD4CGeep0BfCZZ7u8ZwZSLlC7wQSIHmVUq2b+xp46RdSnD59XogTRhg3mzy1MDBnC92DQIOvnfvQRn9ujR+B6H3/M9QYOtH6NsWONB2K7dkRxcWRqFJeZSVShAtddsMD6daPNmjXGb/HIkcB1jx8nSkzk+lWrEhUvztsvv0y0Ywc/rPU9u+QSoqVLo/MZfDF+PPejf3/n+lCuHPfh33+d60NOjvE95Z0xRxI7AsQxLywiygIwHMAsAGsBfEtEq5VSzyml9BKiWQAOKqXWAJgH4GEiOhjKdQPl2fZHYiIHBbz2Wl438PjjofSgYJKZybGvgNy5FszSuzfr4H/9NXCsse3buaxZ0/o1br3VWL38++/sNvzMM0D37oHPS0gA/vMf3v7gA+vXjTaTJnHZvz9QunTguuedx15t2skgI4Pv06OPsn49Lc0IKPrPP7xuJ5hNJVI4af/QuMGV9/Bh/g5Klgz+/TqNo+tAiOgnIqpLRLWJ6EXPvqeJaLpnm4joQSK6iIgaE9HXoV5Tx1kyu8ZAU7QorwVJSuIYUG7LnxxJiIBHHuEAh40b5198aYaKFdkb68wZ9oDyhxYgdrxPEhOBhQvZc27ePGDlSvZAMsPQoWyMnzSJP6eb+eUXLnv2NFe/XDl2Ipg3j722PvvMMKgDLGR//53dWP/910jkFG3cIEC067h+TjiB9sCqXj339+RGXLzGMTJYibOUl6pVgTvu4O3PPgtfn9zEiRP5k0O98grw9tssRMeMsf+j1m6ygcLEhDIDAXhh6BVXAO3b80JRs9SsyQtHz541clG4kS1bOANksWJAhw7mz6tene9Jixb5v7+iRVm4a2GrveyijRsEiBtmIAXFhReIQQGiRxZ2I1wOHMjl+PHs0lqY+PJLVnWULs3urfv28aj1hRf4+LffAh072m9/4EBuf+5cXpORFyLOYgg48xC55houZ82K/rXNMnky36cbbgic59sO/foBZctyOJFoB6ckh114NW6YgaxezaWT98EsMSdAQg2R3Lw5+84fO8brDgoLGRmspiJi28Fbb7HaqWNHQ29+7bWhXaNMGW4X4EyGeQXwrl3A/v38ELM7AwmFLl24nD3bvSvT9SLY1q3D33ZSEjBkCG9Hexaydy/PfMuWZZWbU7hhBqLV482aOdcHs8ScANEjCzsqLM3NN3Opw04XBt5/nxeWNW/Oca7q1+ckURUq8CI7K8EQAzFoEKtLDh3iB7a3usz7j+OE7rd2bX4dOcLGZTeiF7kFyhUfCsOGcfn11/wdBSIzk20nc+YAf/0VWsBRN8w+AHfMQJYt41IEiAsJR5KWzp25nD/ffwiJggQRr1QGOEtjx46sStq/n0eGb70VvnzjSgFffcUPwL//ZmFy6hQfW7iQy2jmgciLzg/jVjWWDrNhxYvQChdeyPfg9Ongg4YHHuCouZ06cVj/22+3L0TcIkCcnoEcOcKDhGLFgAYNnOmDFWJKgGRnG2kidex/O1Styl/uiRMcEqKgs3Yth/xITga6dYv89apWBcaO5e3ly9n2kJnJ9haAH0pOodVY2tPJTRw/zg+2+PjIGlh1SoNRo/xHKj58mO2AABvzk5KAceOAN96wd00tGJ0WIE6HdF++nMuLL84fKdyNxJQA2buX9fvJyaEHl+vUicvCoMbSn+Hqq6P3o+3YEXjME8B/7lygUiVWG5UsyRF8naJDB55tLVnCo0E3sXQpj/Aj/XC55hpet3PyJMeC+/LL/HXGj+eZY+fO/P198QXvf+QRY72QFdwyA6lQgYXhwYP2v/+cnMBrnQKho34XBPUVEGMCJBQX3rxoNdbs2aG35TRagOjPFC1eeolHrMWLG3/W++7jxW9OUaoUh4fPzjZynrgF7Rml0/hGCqWAp54y3t96K/DRR0Yk45wcY8Glnq1cdx3w3HMs4AYMAN59lwNbaqN/MNwiQOLijCCea9YEruuLbdtYwJcpw4FDs7Otna//i6F4O0YVq0vX3f4KFMpEB6Dr3t3/cn6zHD3KMbXi4pwNexAqmZkcEgMgSk93pg+rVhHdfDPRqFEcq8xpXniB78fQofbO37eP72VmZnj7ddttFNWgj6mpRmgPgKh3bw68OGkSv69ePff3lZNDdOONRn0dPmX7dqKDB/1f5/Rp/i8pxQEwneaWW7jvY8daO+/ECaKaNXN//r59ibZuNRdb6+BB45kSLOhoJEBBCmXiBDqMiQ7RHgqlSvHCs5wc4JtvQm/PKb79llUVjRqFZ2Zmh4YNWU3yzDPhM9aHgrch3apR+K23WA1SrRrbyU6fDl+/tAdWpAzoeWnRghctDhnCI/Pvv+cFiTfcwMfvvDP396UU20FeeIFDw1SrxrP+mjVZbfzcc76vs3IlO6PUr++O0PI6zYOeFZnl4Yd5Iewll/C9KlkSmDKFnzfdugX/LX35Jd+HTp2cdWW2hFWJ4/ZXoBnIXXfxqODNNy2I5QB88w2317JleNpzgtat+TP8739O98Q9ZGcbI+9168yfl5HBI26AqFgxLhs1IrrpJg7IGSq67S1bQm/LKvPn82epUoVfLVsS7d8f+Jxp04hSUri+Dmw5ZUr+ejrQ5i23RKbvVvn8c7Ic1HHrVj5HKaI//uB9f/9tpAoAODKyvxn2iRNElSr5v0fRAAUpGm+kXoEEiA4p/v33Fu5qADIyOHeBDpVd0Ni7l//YRYu6Q3XgJnQI/7ffNn/Ogw/yOQ0aGJFl9WvcuND6k5FB50KNu0HNZ5U33uD+N2yYP2z8Qw/xsRdfdKZveZk3j/vTpo35c779ls/xlSfnww+N30Hlyr7D72u1aatWzoXVtyNAYkqFpbMK5k15apekJA75AfgOzeF2RoxgFVy7du5QHbiJ9u251F4xwThxgo3GAHsoDRwIrFhhhIF54onQXEO1+jUlxR1qPqsMH86edqtXc1BSbyK9ONIqOgvgzp2B63mjF8G2aJH/2J13GmruPXtYxXn0qHE8NZXjzQFcuj2AojcxI0CIjD9wONNE3nQTl7Nn88O4oLB/P0cVLlKEAyUKudGeTmZjQi1axOs0WrbkhZBKcdTiRx7hdS179hhhQuwQbftHuElI4ECcAHvfaY8uIPKLI62iI3Xv3m3ei0qnTL7kEt/H+/UzFqeePMleWmvW8Ar+K6/kAUj//tYCZLqBmBEgR49yTKcSJcIbhK5uXRZI+/fziLOg8M47/Ofo2pWN2EJuGjfm1cAbN/KiuWBoQZN3DUt8PIeIL16cFycuXWqvP/oh65ZRuh369ePf2t69bFwGOOaYNla75bMVK8azhOxsFvxm0J+hbl3/dbp04YGGpmFDXsF/6hQb2T//3H6fnSJmBIh3DKxwThGVMhYVFpQ1IXv3Am++ydtPPOFsX9xKfDzQpAlvm4mLFWiNRnKykQagb1/gwAHr/XGbmscOShnrRp59lkfdq1ezp9oFF/Co3C3olf5m1FhZWeZniG3a8ECzbl1W6VWqxB6Q779fMID291EAACAASURBVFae5yVmBEg4FxHmRS/AC3VVOhGHULcydbbD00/zNLpXL/5BC77RMbmCqbGIWBUB+F/k98orPNrcscMQ3lYo6CoszcCB7Ca7bh3bjHQooEgvjrSKVnPr3ByB2L6dhUi1ajzTDEbjxsD69Ty72bOH3ZhTUkLqrmPEnAAJJYiiP666issFC0Lz+3/0UQ6hXrUqT28jERH2iy84DlVcHOuiBf+YtYPs2sUx1sqU8f+AT0w04kR98gnH/rJCYVBhAaxCfvll3h4zxkgu5jbdvxVDultW0TtBzAiQcC4izEuFCmw8O3WKF17Z4bffgP/7P95OSOARSsuWwOjR1to5e5ZVVL6iBK9dCzz0EG+PGSO2j2CYFSDff8/lZZcFVo+2acO/k/37A2dlzEtWlpGpMRK/32jTsyc/oHfs4HtbqpThjOIWrAgQHa4lkP2jsBIzAiTSOuRQ7SD33ssPivvuYx25Xg09apSRHyAYx4+zeqBSJV79u2EDq8N27+a8DW3bsorsqqsMXbTgn7p1+eG2a1fgB4kOJDh4cOD2lAIef5y3H3jAvIvwjh2GiiQpydw5bqZoUbaBlCrFs+0JE9znRm7FBiIzkBgg0gIkFDvIihX8KlsWePVV/jPNnMleKydPchImb7dHXxCxENJrXfbsYWFStSq/rrySI4z26AHMmMEqLCEwcXFGaHl/32tmJq8BUMpcKPwbb+Q1EWfPGjnIg1FY1Ffe3HYbe0ampxuphN2EzEDMEROPEaLI/wnbtmXV09Kl1nOl69S4N97ILoSa8eN5JLRyJXDXXYHbSE3l+gkJwHvvsS6+cmXjVaUKtzF1auEYxUaLYPlBdBynevXMu4ePGsXf86xZxsAmEIXBA6ugYcWILjMQh1BKdVVKrVdKbVJKPRag3nVKKVJK+VjnGZz9+/mhXqpUaImkAlG8OK8BIDISI5nlxx+57N8/9/6kJPaXj4/nIHWB1pnonA13382vjRsN9dXu3ayGKaiugk6iBcivv/r2jNP2EV8rkP1RrhwPFgAjBHogCosHVkGicmVeZLtvX2CHh8xMtk/FxYUvwkVBwjEBopQqAuA9AN0AXARggFLqIh/1SgK4D4Dt3H+rV3PZsGFkwwTYSTKVmckujXFxvl0ZmzfnmQORoT/Py9mznMMa4NwNQvioU4ddLA8dYjtSXuzm6Hj4YfbM+vxzIzujPwqjCsvtFClieGymp/uvt2ULR6BISQk9SV1BxMkZSCsAm4hoCxGdAfA1gN4+6j0P4FUAth1ktQC5KJ94Ci92kkytXcsj2zp1/PuQjxzJSZZ++sn3Q2zRIp5l1asHNG1qvd+Cf5QyhPIzz7AxW0NkrCy2KkAaNTIEx/PP+/aa04gKyxnMGNJ1yt9YnH0AQQSIUqq3Uuoer/dLlFJbPK/rQ7x2VQDeX026Z5/39ZsBqE5EPwbp51ClVJpSKm3//v35juv1FHplcaRo1owN4Vu2GG7DwdAL0C6+2H+dChUM99tHH82v8tAznu7dC1YgtoLCPffwd7BgATs5zJnDhvNff2UDaoUKPFO0ys03s9DftYtjZvmCSASIU5gxpEdygXJBINgM5BEA073eFwPQEkB7AEHMuqGhlIoD8CaA/warS0RjiagFEbVITk7OdzxaqUCLFDFSUZpVY+l6wRZSPfggh8RYvBj4+efcx/SMJ9opaWOFihXZ7RQAnnySVZXNmxuu1iNG2FNfxMXxauz4eF6X4yuF7p49vL6oXDl3hfqIBcwY0r1DJMUiwQRIAhF5y9+FRHSQiHYAKBHitXcB8I6LW82zT1MSQCMA85VS2wBcBmC6VUP6wYOsJoqP9x8pM5xYXQ+yYEHu8/xRsiQLEYC9rTSHD/MMKz6ew7ILkWHQIH517Ai0bm3sv/pq4DG/7h/BaduWVZQAr9DOG9F53TouY9HDx2mszEAiEeGiIBBMgJT1fkNEw73e5h/qWyMVQB2lVC2lVAKA/vCa7RDRUSIqT0QpRJQCYDGAXkRkKcDH7NmsBmjblo2WkUZ77fzwQ/D8D3v38qtkSXMeNjfdxCqqadMMtca8efzQadOGw0QIkSEpiQX3nDkcbWDNGk4HPG1a6Pk5hg0DSpfmtidOzH1Mq1+bNQvtGoJ1tAAJZESXGUhgliil7si7Uyl1J4C/QrkwEWUBGA5gFoC1AL4lotVKqeeUUr1Cadsb7Z1kZpFXOLjgAuC66zgm1gcfBK67ciWXjRubs13UqMHB6M6eZVUKYKjAgs1ghPDSoAHnBg/HoKRiReC113h79Ojci0ajpX4V8qOFQqCBYKzPQALnuwUqAPgDwDwAb3he8wH8CaCi1fSH0Xh5p7TdscNI2bpnj/UUj3bRKTGrVyfKzPRf78UXud5dd5lve/t2ooQEzr2cnk5Uuza3sXhxyN0WHOTECSPn+ejRvC87m6h8ebKcm10ID7t3871PTvZfp2JFrpOeHr1+RQqEO6UtEe0jojZgV9ptntdzRNSaiPZGSKaFjY8/ZvXOdddxfKho0a4dUL8+6051mlNfaKOplUikNWpwMDoiNuxu3szqDzteQIJ7KFECePdd3n7/ff7dpqZyXLSaNWMzTIbTVKjAjjH79/teTHj2LC80jIvjWWQsEsyNN1EpdT+AvgDOAPiAiOZGpWchkpFh+NkHCwMSbuLigBdf5O3nn+fEOXnJyAAWLmTVldVQ1kOHcqmFU6dOBTNPtpCbnj1ZWGzZwtsffsj7+/QR92wn8F5MqG0d3vz7Lw/kKlaM3f9fMBvIZwBaAFgJXjFuMbi4czz6KBuoL77YGe+ka6/lhER79/ICtLxrNxYtYl1306ZA+fLW2u7cmWdVAD9Y/K0hEAoWRYoA99/P2z/9xEb7IkWiPwASDALZQWLe/gEgmNy8iIgaA4BS6n8I0XAeLY4eBT79lLfHj3dm9KYUZ6Hr2JEz0B0/bhi+iYw8H9pry2rbX33FD5nzzzcy5wkFn3vv5Znpd9/x+zFjeLGh4AxagPjyxIr1RYRAcAFyLsACEWWpAjKPnjmTVUTt2jkb2qNDB1478NlnrG7Kaw9JSuLQ3naIjwd6+wr8IhRo4uL499KiBbt265mm4AzVqnHpawaihYquE4sEEyCXKKWOAdCSI8nrPRGRyQDW0WXWLC579HC2HwDPgKpX54eCN6VLc1iMWB69CL4pUSK0xYlC+AikwtILDKtXz38sVggoQIioSLQ6Ek60ANGhJpzm+ef5JQhCwSKQCksESBABopRKBDAMwIUAVgD4lHgBoGs5dYo9JipWDBygUBAEIRiBVFgiQKx5YXUHLyR0NQcOcNmzp7g+CoIQGqLCCkyh88LS6WSHDXO2H4IgFHy8BUhODjs5AJwXRq8NiWU33mAzkFxeWBHuS1jIzua4QbIyWxCEUElM5FD6WVm8Il2zZw8LlIoVOb99rBJMgFyilDrmeR0HcLHe9nhjuZK773a6B4IgFBZ8GdJFfcUEi4VVhIhKeV4liaio17YrXXiLFAFuvNHpXgiCUFjwZUjXwkQESCGjUSNeoCcIghAOfBnSZQbCFDoBEqtBzQRBiAyiwvJPoRMggiAI4cSXCksLkFgOYwKIABEEQQiIzED8IwJEEAQhACkpXG7ebOzbsYPLWBcgYjEQBEEIQO3aHP162zbg5EnORLh3L68RifVgqDIDEQRBCEB8vJFSeN06YM0a3m7QgJcNxDIiQARBEIJw0UVcrl5tCBC9L5YRASIIghCEhg25XLMG+Ptv3pZo3w4LEKVUV6XUeqXUJqVUvhQ6SqkHlVJrlFIrlFJzlFI1neinIAixjfcMJDWVt1u2dK4/bsExAaKUKgLgPQDdAFwEYIBSKu+kcBmAFkR0MYDJAF6Lbi8FQRCMGcj8+cDSpZwqQgK2OjsDaQVgExFtIaIzAL4GkCvLNxHNI6IMz9vFAGJ82Y4gCE7QoAFQvz5w4gRAxDnrS7kyGmB0cVKAVAWw0+t9umefP24H8HNEeyQIguADpYCnnzbed+vmXF/cRIFYB6KUugWcGfFKP8eHAhgKADVq1IhizwRBiBUGDACKF2dDuqSMYJwUILsAeK/jrObZlwulVCcAIwFcSUSZvhoiorEAxgJAixYtKPxdFQRBAHr35pfAOKnCSgVQRylVSymVAKA/gOneFZRSTQF8BKAXEe1zoI+CIAjuICeHXcD27HG6J+dwTIB4UuQOBzALwFoA3xLRaqXUc0qpXp5qrwM4D8AkpdRypdR0P80JgiAUbiZOBFq1AurVA44edbo3ABy2gRDRTwB+yrPvaa/tTlHvlCAIghuZM4fL48eBuXOBa691tj+QleiCIAgFA70EHgBmz3auH16IABEEQQgXqanAe+/lzj4VDrKzjSBcALB4cXjbt4kIEEEQhHDwyy/AZZcBw4cDl1wCbN0avrb37weysjiGfFwcsGIFcOpU+Nq3iQgQQRCEUCECnniCPaWUAg4eBEaMCF/7u3dzWbcux1XJygKWLQtf+zaJbQFy+jTwzz/AmTNO90QQhILMunUcJKtcOc48dd55wI8/8vMlHGjX3cqV2RMLAP76Kzxth0DsCpCjR3ma2aQJ0LgxsH270z0S/DFlCgcfmjDB6Z4Igm+0Ufvqq4EaNYDbbuP3H30Unvb1DKRKFeDSS3l70aLwtB0CsStAXnsN2LCBtzdsAK68EjhwwNk+Cfk5dAjo359Hd7fcAqxc6XSPhMLCkSMc4OrZZ3k7FBYs4LJDBy6HDuXyyy85AmOoeM9ArrqKt3/+GcjI8H9OFIhNAUIEfP01b0+fzoH9t28HBg3iY4J7mDyZk1BrHnnEub4IhYfMTKBzZ+D554FRo4AePdiuYJcVK7jUMd4bNQLatOE1G/pZEwr7PIE4KlYELriAn1knTwI//RT4vAgTmwJkwwZgyxYgORno3p0fUuefz1/GpElO907wZt48Ll94geNnz5wJPPBAeEZ1QuzyxhtAWhpQoQJQtizwxx/AN9/Ya+vUKWDjRk6Q7p3n9s47ufzkk9D7e/Agl+XLc3njjVza7XOYiE0BsmQJl1dcwV96jRrASy/xvldecW4WQgR89RWH/Xz44fD7khdEtL97r17Ahx+yh8uYMUCfPjJbFOxx4AD/zwEOD/L667z97rv22lu/nr2v6tYFihUz9l9/PRvTlywx1OWh9BkwBMgNN3D500+OqrFiU4CkpXHpnZNy0CD2oFi2zMhZGW2efRa46Sae8o4ezQb+1aud6YsbOHiQPVpKlOCR3YABPFIsX57DOrhkNW4+UlOBV19l9ajY1dzHCy+waqlrV7Yn9O8PJCXxYCU93Xp7mzZxWadO7v3Fi7MQAdgWEgp5BUiNGuyNlZHBs3KHiE0BovWVTZsa+xITw+85YYU1a1gfGxfHo6OOHflH07kzLyKKRdat47J+fZ4pArxQ6667eHvqVGf6FYiZM4HWrYHHHuO438nJnH3IRRFUY5p//wXef59nsnoWUqIE0KULb//yi/U2N2/msnbt/MduvZXLL77gWYpd8goQwBBODqrdY1OA6BFD3bq59//nP1xOnsxrRKLJmDH8Axs6FHj0UeCHH4DLL+cHz113uVtdk5HBo+5w2yXWr+eyfv3c+3v25NLOnz2S5OSwfSY7mx9IV17JA5OZM7nP3s4AgjN8+il/D716sRu/RntP/f679TYDCZD27YHq1Xkm/ccf1tsG+L+vBUi5csb+vn25nDmTf3MOEHsCJCODbQvx8TwN9KZePaBZM+DYMXaRixZZWcB33/G2Xr2alMTT3vPO42Nue1hqtm1jQdyqFbsYDh/O6qVwPCz1DKRevdz7mzRhXfPmzaG7X4aTBQu4z9Wr8wBg/ny+PzVrshvyyy873cPYJjsbGDuWt/UsVtO2LZd21lZs28ZlrVr5j8XFGfaKGTOstw3wM+v0aR6MFC9u7K9dm19Hjhhq+SgTewJkyxYua9UCivqIZt+/P5fRXLS2eDGvd6hbN7cXR0oK8NRTvP3YY6FNgSPFiBEskEuU4BnIe+8BnTqxfSlU1ZuegeQVIPHxxujRBeEczqEHATfdxH0E2O1y3Djefvllezp2ITykpbG7fo0arBr2plEjICGBtRPHj1trVy/yq1rV9/Hu3bm0OyjV/6PkZFa9eaM/x6+/2ms7RGJPgGj11YUX+j4+YADr26dNAz7/HJg1K/KqBz16uNJHyvcRI/iHuXy5/RFMpNi4kUfaSUksmFesAB58EKhWjUM4DBsWWvv+VFgAzxSB3CGunWbhQi67dcu9v0MH1lefPs1rDsLBG2/wTKdNG8Or0CmIeE2Cm9WsAPDbb1xefTXPDLxJSOAYU4BhIzWL9ypxX1xxBQ+wVq60N4DQa0AqVMh/rJMnZZJDDiUiQPJSrRowcCBPdwcNYk+NXr0i++dYvpxLb6O+JikJ+O9/efuddyLXBzto49311/OPu3FjfrD9+Sf3e8oU639GzdmzrKJSKr93C2As2Fq61F774ebUKf6scXFG37x56SUemIwblzsstx0mTQIeeogfRn/+ybayr74KrU27zJrFC9vOO48FvR0bQrSYP5/L9u19H9f/Pyuz2lOngMOHecbpbZ/wplgxY/X4rFnm29YEEiAdO/Jv7o8/HFkbJQLEF2+9xaqsBg34Tz9zJs9IIoX+wTZp4vv4bbfxCGbOnNAfPuFk8mQutY5XU60aMHgwb3/+ub22t21j21D16iyM8qJnIG4RIMuW8aCjYUN+mOalTh12kMjJAR5/3P51cnKAkSN5+9lngfvu4+vecgvw/ff227XD5s3sabZtGz/ENmxglYr+XbiJrCwj3IivmT5g/P/0gM4M3iFG8s5qvNEzhblzzbetCSRAypblOHFnzzoivEWA+KJkSR7RrVkD/N//8T5dhpvMTL6OUsDFF/uuU6YMPyAA9iJxA5s380OzZMn8+mTA6O9XX9nzENGGyQsu8H1c66w3bGCnB6scOQLccQc/2AcMCD13g147pCOl+uLpp3kgMH26oe6yyu+/s+qwWjUOHz5mDNvJcnJ40PPnn/batcNLL/Hvt18/HokPH86Rrfv140HP11875h2UjxUreIReu7Z/W4UWIFZmIMHUVxrt5TVvnnVthrcNxBcO2kFEgARj8GB+SC5YEJlAfqtX8+ioXj1+uPhDr1H58kt3uINOnMhlr17sHZKX1q3ZCWD3bmPkZwUtQFJSfB9PSGCVGWDdkJ6VBVx3HYeY2LSJH3SNG4dmY9KhtQMJkEqV2M0X4IWGdtCj+1tvNZxAnn2WheHp0+wuHOqqZzOcOcMqSoDtOgkJwNtvs1AhAsaPZ8F8yy2hr384eTL0/mpVaosW/utox4xVq8z/x8wKkIYNWQDs2cMDACsEmoEAIkCiRmYmsGMHq6Vq1jR3TsmSxmKgF14Ivy0kmPpK06oV65j37rWnRw0nRIaXmp5p5EUpHokCPOK2ig6vH+h70g8Dqy6Mr77KqoQKFdg9+rrr+CHVp499TxkzAgRgp4hixdj5wM6DXgfP69XL2KcUL47r0YNX73ftygvmIsncuTyLa9iQVb26H48/ziqgl17i/87XX7MNz+r/ZsUKNj4nJ3Ocup49WUVnN+ChjuigDeW+KFWKvTPPnDEcOIJhVoAoZdhedHw3swQTIK1bs9p09Wp+vkWR2BIgW7fyD7lmTR4xmeXBB3mU/e23rOt8/fXwCZJABnRvlDLsCuPH+65DFJ0FkKmp/AerUMHQ7fpCL/ibMcP6/Qo2AwHsGdIPHjRG/xMm8Oht0iRevJmTwyN5q7GFDh3imUxSUuAHFMD3TAvdt96ydp30dP4Nly6dOwwPwLORb77h/Vu3sndWJO1l/uxfAI/kH3+cIwUULcpqtiFDzI/qFy3inBeLFvE9PXuWBW6fPiyg7TwkV63iMtj3owdyZhNBmRUgQG41lhW0APGnwkpIYM8yIPqemkTk2AtAVwDrAWwC8JiP48UAfOM5vgRASrA2mzdvTn6ZMYMIIOrSxX8df0ybRlS5Mp8PEE2YYL0NX1x+Obc3a1bwuunpREoRJSQQHT6c+9ihQ0Rt23JbvXoRnToVnv7lJSeHqHt3vs5DDwWue/Ys0fnnc921a61dR9+XefP811m6lOvUrWu+3cce43Ouvjr3/uxsombN+Nhzz1nr68yZfF6bNubqr1rF9YsXJzp40Px1vvuOz+vc2X+dffuIWrTgeqVKEc2da759s5w9S1S+PF9j1arAdWfM4M8JEN13X/C2z5zh7xMguuUWomPHiP79l2j0aKKaNXl/9epE69ZZ63P16nzu+vWB640axfUefthcuzffzPXHjw9ed+1arluxIv+PzNK0KZ+Xmuq/zmefERUtSvT44+bbzQOANLL6DLd6QrheAIoA2AzgAgAJAP4BcFGeOncD+NCz3R/AN8HaDShA/u//+CPffbed+8tfum6jRg3+I4VCdjbReedxe/v2mTunY0eu/7//5d5/3XWGcAOInnjCen9+/JFo+HCi33/3X+fzz7n9kiX5jx2MW27h+q+/bq0vVavyeVu2+K+TmUmUlMT1zPRl716iEiW4/pIl+Y/PncvHypQhOn7cfF+fe47Pu/9+8+d06cLnvPKK+XMefZTPGTkycL2TJ4muv57rnn8+0a5d5q9hhnnzDMFt5kG4aBFRXBy/tm0LXHfqVG67dm3+fr05fJiFNEBUujRR7948mLnjDt/fp+bIET6nWDGirCxz1zc7yOzQgev/8kvwujk5RJUqcf3Vq821T2T8F7Zv91/n5En+nCFQ0ARIawCzvN4/DuDxPHVmAWjt2S4K4AAAFajdgALknnv4I7/5po3b6yE72xghff21/XaIiDZs4HaqVjV/zscf8zlXXWXs0yPx4sWJJk40ts08VDVz5vDsBuCRzOLF+evs2MF/XIDo00/Ntfv111y/XTvzfcnM5L7ExfGINBD6QTxxYvB277+f615zjf86+gH1zjvm+3vNNeb7oPn5Z+O7D/YZNe3b8znTpwevm53NsyyAqGtXayPeYAwfzu0++qj5c266ic958snA9W64geu99prv4ydO8OfxHiwB/Ht57z3f5/zxB9dp0iR4P7ds4boVKpi7Z/Xrk6mZmKZ/f67/7rvm6ufkEMXH8zkZGebOsUlBEyDXA/jE6/2tAN7NU2cVgGpe7zcDKB+o3YACRP+hzPwBA/HBB9xOy5ah/TG/+Ybb6dHD/DmHDvEPKi6OaPdu3nfXXdzOvffy+169+P3w4ebazMkhuvRSPichgcuaNflamqwsok6d+FjPnuY/95EjLJCKFDGvrtm0ic7N8oLx2mtcd8iQwPW8R8F//+2/3uTJdG4EHGy0SsT3oUIFPmfTpuD1vc9r0MC84MnKMmare/aYu8aePTybAlgFq8nIIBo2jO9vt278MHv5Zf69BBtJHzlitBlIpZKX2bP5nAYN/NfJyjLaDnQvs7OJ5s8nmjKF/8t6YFCkCNGKFfnr60HXzTcH72dOjtGH9PTg9UuV4rpmf9sffcT1r7vOXP3Dh7n+eeeZqx8CMStAAAwFkAYgrUagh06tWvyR16yxdYPPcfKkodtfuNB+O48/TqZGZXnRAmLMGH4Y6FnBP//w8ZUr+UFZtCjPcoKxcCGdU3ccPsyCEWAVQVYW/6keeYT3lStn/gGm0Wq3L780V18/bNq2DV532TI6pxf3J9SOHye64AKu99hjgdvLyjJ+J1OnBr/+9u3GvbM6mNAPk0svDV5340ayPFslInr7bT6vVi3+reTkEA0cSPlG8N4vf7PLEyeM2ZaVGSURz7L0g9nfbzItzeirVYYN43NvuCH/sfvu42MvvWSurauuyi90fXH8OJ1TjZn97rXWoVw5FoTBWLPGGNBEmIImQKKrwsrI4Glu0aL5dat2GDmSb9+119pvQ0/FJ0+2dp5WC116KT+UATacenP77by/Q4fgI2ktkLTdZMsWQyi1bEl0xRV0boQ3Z461vhIRvfEGn3/rrebqf/KJ+frZ2YZB15+B9M47+fgll5j77seMMS/AJk2ic2oiq5w8aX40//33XC+QAd0XZ88SNWrE5z77LM82tIrzp5+IPvyQBcr997NtEGA70ebNudvZvp1VQNr+YNWITWSosd54w/dx3beBA623nZ5uDJry2hP1zNms5kHbmoIN7LQwsCLwcnIMm4Ye8AVCqzq9VdYRoqAJkKIAtgCo5WVEb5inzj15jOjfBmvXrwDRI9X69e3e39zs2WOoe6xM5TU5OeyNAeT/swbj5EnDGKx/jB98kLvOv/8a7T/9tP+2li/nOomJuW0m8+fzKEmPSkuWtG/zWbmS2zDrffLUU1z/qafMtX/jjeRXr7x2LQ8c4uO5H2Y4dsy4v8G+m4cfDn6PA/Hf//L5gwcHrvfKK1zPjCdTXrTR2/vlT22mdfTt2vHAY/9+nilr9dmFF9qfweuBz5VX+j6uZ0b+bBnB0AOyvA4m2nBt9n9mdlAwfz7Xu/xya/3UjiVvvx28rlaX3367tWvYoEAJEO4vugPY4FFNjfTsew5AL892IoBJHjfevwBcEKxNvwJkwgT+uH372r7B+dAPj8svt66+2LyZzy1f3p4dRY/m9Ijw6NH8dbRhXCnfbsLZ2YYR2pcH0dGjrGeeOjW/27AVcnKIqlQh06OuW2/1/SDwx//+x/W7dct/TP9Z77zTWp/1/Q2m9rjySq73ww/W2tds2sTfT7Fi/LD2h364fvihveu8+CKP0IsUIXr1Vf/1DhwwbDrNmhkuuADPVL3tYlY5csSw3x04kP+4Nkj/9Ze99t9/n8/3ti/s30/nZlVmVEZERFu3mvtvaocVX2qzQOgZtplnkZ4NPf+8tWvYoMAJkEi8/AqQJ5/kjxvMBdIKR44QJSdzu1ZH59odtlcve9f+80/DO2PMGP/1tItp+fK5jYJnz/KoV+vvrdo1rKKvZcadV6vMZs821/bevfxQSkjI7cq4+SC6uwAAEV5JREFUcaOh1ti61Vp/tcro4ov91zl92nAjNuuG7Qu9riaQS6+2SwVysQ7Gvn2BhZRm7lyekWrB0bUrezKFA61OymsP0662CQl8X+2gH/ylShmebXr21aqV+XZycozZdyC349GjydasUDuJmLGD9OvHdT//3No1bGBHgMTOSvS1a7n0TtgUKqVLAy++yNuPPMIB5cyiw5H4iwwajMsu49AFCxZwRFZ/jBzJ6VUPHODVquPGAT/+yCuwx4/nDGfffcdxmiKJlZzTgVKE+qJCBc4od+aMEeoD4HAaOTkclj/QinZfXH01h7ZYscLoT17++ou/cx3nyC7Dh3P50Uf8yM4LUXh+v8nJuXNq+6NDB4408PnnvIL75585XEY4uOYaLn/4Ifd+HY5GZ5u0Q0oKh1U5dszILKhXoOu4aWZQygiTEyjKgZVV6N5ccAEHwzx40Aix4g+dlTNv+m2XEDsCRId10HF7wsWQIRxFd8cOzoVhhpMnjT+Qd0wjq9Spw/GCAhEXxwEYL7yQf6xDhvCfeP58zpb3yy/+8yOEEx3y5PffAwfHy8jggHPx8RzK3Sw6P/QXX3C5ahU/AIsUsRc+vVgxIzGUv/AQOny6zvVgly5d+IGydavvlKo7d3Ik2eRk/zknwk2NGhwDLljoD6toATJzZu64VjopVrBYYsHQ35keqGgB0qiRtXZ0mJxAcdbsChDvuFg6R4kvsrMDJ1VzAbEhQI4f5y+iaNHwfxFFinCsH4BTlu7dG/ycd98Fjh7leD9mowKHQnIyj6TeeYcDB151FSck+ucfTkYUDZKTedaUmRk4uKJOOZySwvfWLDfdxHGTfv4ZeP55jpuUnQ3ceaf5mUxetHD3JUCOHTOElU6DbJciRYz4WL7yp5iN41QQqF2b/4NHjnASJM3ixVyGOtPJO9PVEbStChAzgTrtChDAECBz5vivs307/1+qVGFthxuxqvNy+8unDUSvK2jZ0rpi0CzaP37UqMD1Fi82bBc//hi5/rgR7abZvbv/OtOmGXp3qzz/vKG3BziG0LFj9vt76BAbnYsWzb1QLD3diDd16aXhWeWt/f1Ll85vA3j1VT42YkTo13ED2vNMx5vKyTFcsa16JObl5El2SFCKF9pqb7q9e621s2MHn1e2rP/v98ILuY7VOG9ERDt30jl3an8rzLUdrlMn6+3bAGID8YPOAtamTeSu8eCDXH7wAevifXHkCIc4P3uW7Rbdu0euP26kXz+eBc6aBeza5buOVfuHN088wbOsa6/l0PsLFnBIcbuULcspQ7OyjPS9hw/zKDctjWdJn33GKolQadCA9f9Hj+YP1293FO1W8tpBtmxhG12FChxOPRSKF2d7GBHwyiusLq1Tx38odH9Uq8az5sOHfUf/JTKyEdqZgVSrxrOcjAz/+cz1927FfhNlCr8AIWIjMWCEF48E7dvzH3zvXuNhk5dHH+UfY6tWwGuvRa4vbiU5mR/u2dmG80FeQhEgcXFskJ4yhZ0HAiXoMsvAgVy+9RY/jHr3ZnvaRRdxWPt69UK/hkarwr7+Ovd+O4ZgN3P55aySWbuWhYdWZV12WXiEsVY9vv02l+3aWW9DKcNe6itvy/Hj/HsoUcL+IKV3by79pcsWAeIwx46xYXv9evYysvNDMotSwL338vZrr+VP5fnbb8DYsWwcHjfOWj6SwsSoUfyg//hj35nZtACJhm3IDP36cf6YtWs5F8yCBZwSdeZMcx5NVrjxRi6//95wNMjKMjywCoMNBOD/gM5f8eOPoXsk5mXAgNyeXP6SngWjTh0ufQkQb/uHXaHXpw+XM2b4Tv0rAsQhiDhPdHIy8PDDvO+ZZ/iHG0luuYW9V1asYJdMzeHDRkraJ54IrytxQeOiizgxVlYWcNdd+d1Wtdui/vM6TUICq6mSknjUWbEiG+qteIiZJSWFHSsyMvjBCrCQzcxkIVaqVPiv6RRajTVpkuF6HS6VbvnyPPuIjwduv92+YNKus74GOlqAVK5sr22ABwQXXMA5z70dCgD+ztev58GWm58XVo0mbn81b96cV+tqQ+pll7HxNpzhrAOhk/4kJnJI8MWLjWCCzZvbXyRVmNi/31io9cUXxv5Dh3hfUpK5SLjRZOdOTo4Uyop8M+h8MzrG2tixlG91dWFg/352TvB2eAj3f9TsynN/6NwgviIcfPEFH+vfP7Rr6CCleSMl6NBLVpKlhQhkJTqheaNGhudFuLIGWiEnxwje5/0qXz5wcqRYY9w4vi/JyUZYC53QyUx02sKKzjpZrBjHJgsU56ug88ADxv9jyhSne5MfnTnywgvzH3vpJT4WLCtnMHScuDJlcmcR1ZEqojhwsCNACp8Ka/Nm1h/feCOvDYg2SrEn1oQJPCVv1IjVV6mpoXuYFCYGDWLHg/37eTsnh+1EQP5837FE1ars7JGZyb8bvVixa1dn+xUJXnuN1YMzZrBzhduoXZv/z1u35ves1J5ZNWqEdo1GjYCmTdlD03t1fgGwfwAohDMQgMMr+wrWJriLbdvYzx4geuEFIyf5jBlO98xZVqwwIj0D5tOrCuEnJYW/g7zh63v04P1mcsYEQ6ste/Y09unIwt99F3r7JoHMQMCG8wULohfyQbBPzZocZgUAnnwS+Ptv4PzzORZTLNO4MTBxomFU//hjp3sUu2hnjryG9HDNQAD2GouP5xnIqlU8bPCODeZiCp8AqVGD1QBCwaB7d178FxfH6oLXXw/P+o2CznXXsepk8eLwPKQEe2hPrLyuvDt3chkOb7yKFTnkDhF7aW7YwAsrK1Vyvdq7qNMdEAQMH84PzOxsXqErCG7B1wzk2DG2WSQmhm8t0MiRvD5sxgwgPZ33tW0bnoWVEaTwzUCEgknlyiI8BPehZyB6fRJgzD5q1AjfA75SJSMo67JlXN58c3jajiAiQARBEPyhF/HpaABAeNVX3tx+uxHNondvoEeP8LYfAUSACIIg+KN6dbbJ7d3LCaCA8BrQvVGKY64dO8bxsYq638IgAkQQBMEfcXFGUEWdlE7Ha7Oa5dIsoUSQjjIiQARBEAKhg1hqAaKzBIYzEnMBRQSIIAhCILQdROcvFwFyDhEggiAIgfAWIJmZwKZNbK/QHloxjAgQQRCEQDRtymVqKkdLyMrivO7FizvbLxfgiABRSp2vlPpVKbXRU5b1UaeJUupPpdRqpdQKpdSNTvRVEIQYp2pVXhF+/LgRVubSS53tk0twagbyGIA5RFQHwBzP+7xkABhIRA0BdAUwRilVJop9FARBYDp25HLcOC6vuMK5vrgIpwRIbwCfebY/A9AnbwUi2kBEGz3buwHsA5ActR4KgiBoBg40tuPjgb59neuLi3BKgFQkoj2e7X8BVAxUWSnVCkACgM1+jg9VSqUppdL2798f3p4KgiC0bQsMHcopjt98EyibT+sekygOAx+BhpWaDaCSj0MjAXxGRGW86h4mIp/fiFKqMoD5AAYR0eJg123RogWl6VDIgiAI4eTMGRYihRCl1FIiamHlnIitlSeiTv6OKaX2KqUqE9Eej4DY56deKQA/AhhpRngIgiBElEIqPOzilAprOoBBnu1BAL7PW0EplQBgKoDPiWhyFPsmCIIgmMApAfIKgM5KqY0AOnneQynVQin1iadOPwDtAAxWSi33vNydnksQBCGGiJgNxCnEBiIIgmAdOzYQWYkuCIIg2EIEiCAIgmALESCCIAiCLUSACIIgCLYQASIIgiDYQgSIIAiCYAsRIIIgCIItRIAIgiAIthABIgiCINhCBIggCIJgCxEggiAIgi1EgAiCIAi2EAEiCIIg2EIEiCAIgmALESCCIAiCLQpdPhCl1HEA653uh0soD+CA051wCXIvDOReGMi9MKhHRCWtnBCxnOgOst5qUpTCilIqTe4FI/fCQO6FgdwLA6WU5Ux8osISBEEQbCECRBAEQbBFYRQgY53ugIuQe2Eg98JA7oWB3AsDy/ei0BnRBUEQhOhQGGcggiAIQhQQASIIgiDYosALEKXUNqXUSqXUcu2GppQ6Xyn1q1Jqo6cs63Q/o4Gfe/G6UmqdUmqFUmqqUqqM0/2MBr7uhdex/yqlSClV3qn+RRN/90IpNcLz21itlHrNyT5GCz//kSZKqcV6n1KqldP9jAZKqTJKqcme38BapVRrq8/OAi9APHQgoiZe/tyPAZhDRHUAzPG8jxXy3otfATQioosBbADwuHNdizp57wWUUtUBdAGww7luOUKue6GU6gCgN4BLiKghgNGO9i665P1dvAbgWSJqAuBpz/tY4C0AM4moPoBLAKyFxWdnYREgeekN4DPP9mcA+jjYF0chol+IKMvzdjGAak72xwX8H4BHAMS698hdAF4hokwAIKJ9DvfHSQhAKc92aQC7HexLVFBKlQbQDsD/AICIzhDREVh8dhYGAUIAflFKLVVKDfXsq0hEezzb/wKo6EzXoo6ve+HNEAA/R7lPTpHvXiilegPYRUT/ONu1qOPrd1EXQFul1BKl1G9KqZYO9i+a+LoX9wN4XSm1EzwTi4VZei0A+wGMU0otU0p9opQqAYvPzsIQyuQKItqllKoA4Fel1Drvg0RESqlYGW3muxdE9DsAKKVGAsgCMMHRHkYPX7+LJ8Dqq1jD170oCuB8AJcBaAngW6XUBVT4/fp93YvrATxARN8ppfqBR+WdHO1l5CkKoBmAEUS0RCn1FvKoq8w8Owv8DISIdnnKfQCmAmgFYK9SqjIAeMqYmJ77uRdQSg0GcA2Am2PgAQHA5724Ejzq+kcptQ2syvtbKVXJsU5GCT+/i3QAU4j5C0AOOLBgocbPvRgEYIqnyiTPvsJOOoB0IlrieT8ZLFAsPTsLtABRSpVQSpXU2+DR5SoA08E/CnjK753pYfTwdy+UUl3BOv9eRJThZB+jhZ97kUpEFYgohYhSwH+gZkT0r4NdjTgB/iPTAHTw7K8LIAGFPCptgHuxGzzAAICOADY608Po4fnd71RK1fPsugrAGlh8dhZ0FVZFAFOVUgB/lolENFMplQqekt8OYDuAfg72MVr4uxebABQDT9cBYDERDXOum1HB571wtkuO4e93kQDgU6XUKgBnAAyKgdmpv3txAsBbSqmiAE4D8GU/LIyMADDB81vYAuA28KTC9LNTQpkIgiAItijQKixBEATBOUSACIIgCLYQASIIgiDYQgSIIAiCYAsRIIIgCIItRIAIQhA8UUvv9npfRSk1OULX6qOUejrA8cZKqfGRuLYgWEXceAUhCEqpFAA/EFGjKFzrD/CiT7+L+pRSswEMIaJYiygsuAyZgQhCcF4BUNuTL+J1pVSKZwEelFKDlVLTPLkTtimlhiulHvQEqFuslDrfU6+2UmqmJ4jfAqVU/bwX8awIz9TCQyl1g1JqlVLqH6XU715VZwDoH/mPLQiBEQEiCMF5DMBmTw6Jh30cbwSgLzgo4YsAMoioKYA/AQz01BkLDlzXHMBDAN730c7lAP72ev80gKuJ6BIAvbz2pwFoG8LnEYSwUNBDmQiCG5hHRMcBHFdKHQXPEABgJYCLlVLnAWgDYJInjAbA4WXyUhkcYluzCMB4pdS3MIL9ARzgrkoY+y8IthABIgihk+m1neP1Pgf8H4sDcMST8S4Qp8AJjQAARDRMKXUpgB4AliqlmhPRQQCJnrqC4CiiwhKE4BwHUNLuyUR0DMBWpdQNAKCYS3xUXQvgQv1GKVWbiJYQ0dPgmUl1z6G64CiyguAoIkAEIQieUf8ij0H7dZvN3AzgdqXUPwBWg1OH5uV3AE2Voed6XSm10mOw/wOAzqTYAcCPNvshCGFD3HgFwUV4MsPNIKLZfo4XA/AbOLNelq86ghAtZAYiCO7iJQDFAxyvAeAxER6CG5AZiCAIgmALmYEIgiAIthABIgiCINhCBIggCIJgCxEggiAIgi1EgAiCIAi2+H8vJ9JcUeNiVwAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "fig, ax = plt.subplots()\n",
        "t = np.arange(0, len(ppg_filt))/segment_data.fs\n",
        "\n",
        "ax.plot(t, ppg,\n",
        "        linewidth=2.0,\n",
        "        color = 'blue',\n",
        "        label = \"original PPG\")\n",
        "\n",
        "ax.plot(t, ppg_filt,\n",
        "        linewidth=2.0,\n",
        "        color = 'red',\n",
        "        label = \"filtered PPG\")\n",
        "\n",
        "ax.set(xlim=(0, n_seconds_to_load))\n",
        "plt.xlabel('time (s)')\n",
        "plt.ylabel('PPG')\n",
        "plt.xlim([50, 60])\n",
        "\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rJ4aQtFvpHWI"
      },
      "source": [
        "<div class=\"alert alert-block alert-warning\"><p><b>Note:</b> The PPG signals in MIMIC have already been filtered somewhat by the clinical monitors used to record them.</p></div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iJB6_3eUpHWI"
      },
      "source": [
        "<div class=\"alert alert-block alert-info\"><p><b>Further work:</b> Several different types of filters have been used to filter the PPG signal (e.g. Chebyshev filter, Butterworth filter). Have a look at <a href=\"https://doi.org/10.1038/sdata.2018.76\">this article</a> for examples of several filter types (on pp.8-9). Which type of filter do the authors recommend? Can you re-design the filter above to use this type of filter?</p></div>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "b5UDOyp8pHWI"
      },
      "outputs": [],
      "source": [
        ""
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.8"
    },
    "toc": {
      "base_numbering": 1,
      "nav_menu": {},
      "number_sections": true,
      "sideBar": true,
      "skip_h1_title": true,
      "title_cell": "Table of Contents",
      "title_sidebar": "Contents",
      "toc_cell": false,
      "toc_position": {},
      "toc_section_display": true,
      "toc_window_display": true
    },
    "colab": {
      "name": "signal-filtering.ipynb",
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}